14D-R141  045  AN  APPLICATION  OF  THE  H-FUNCTION  TO  CURVE-FITTING  AND 
DENSITV  ESTIMATIONS)  AIR  FORCE  INST  OF  TECH 
HRIGHT-PATTERSON  AFB  OH  SCHOOL  OF  ENG I.  . 

UNCLASSIFIED  C  D  BODENSCHATZ  ET  AL  DEC  82  F/G  12/1 


fife* 


AFIT/GOR/ OS/ 33D-4 


i-\>>S 


P 

tod 


in 


< 

1 

Q 

< 


a 

*-U 


AN  APPLICATION  OF  THE  H-FUNCTION 
TO  CURVE-FITTING  AND  DENSITY  ESTIMATION 

THESIS 

Carl  D.  Bodenschatz  Ralph  A.  Boedigheiraer 

First  Lieutenant,  USAF  First  Lieutenant,  USAF 

AFIT/GOR/ 0S/83D-4 


DTIC 

ELECTE 

MAY  1  5  1984 


Approved  for  public  release;  distribution  unlimited 


84  05  15  058 


AFIT/GOR/OS/83D-4 


AN  APPLICATION  OF  THE  H-FUNCTION 
TO  CURVE-FITTING  AND  DENSITY  ESTIMATION 


THESIS 


Presented  to  the  Faculty  of  the  School  of  Engineering 
of  the  Air  Force  Institute  of  Technology 
Air  University 

In  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 
Master  of  Science  in  Operations  Research 


Carl  D.  Bodenschatz,  B. A. 
First  Lieutenant,  USAF 


Ralph  A.  Boedigheimer,  B.S. 
First  Lieutenant,  USAF 


December  1983 


Approved  for  public  release;  distribution  unlimited 


Acknowledgments 


We  want  to  express  our  sincere  appreciation  to  our 
thesis  advisor,  LtCol  Ivy  D.  Cook,  Jr.,  for  his  inspiration, 
motivation,  and  guidance  throughout  this  effort.  Our  thanks 
also  go  to  Dr.  Dennis  W.  Quinn,  our  reader,  who  gave  us 
constructive  suggestions.  LtCol  James  N.  Bexfield  deserves 
credit  for  shaping  our  education  in  probability  and  statis¬ 
tics.  A  special  note  of  thanks  also  goes  to  our  typist,  Liz 
Walls,  who  agreed  to  type  this  thesis  in  spite  of  the 
complicated  notation. 

Most  importantly,  we  want  to  thank  our  wives,  Debbie 
and  Becki,  for  their  support  throughout  this  educational 
endeavor.  Without  their  patience  and  encouragement,  this 
experience  would  not  have  been  bearable. 


Accession  For 

NTIS  GRA&I 

DTIC  TAB 

B 

Unannounced 
Justification _ 

□ 

Bv.  .  . 

Dint  ri but  i  on/ 
Availability  Codes 
!Avr- '  I  and/or 
jDist  I  Special 


M 


Table  of  Contents 


Page 


Acknowledgments  .  ii 

List  of  Figures .  v 

List  of  Tables . vii 

Abstract . viii 

I .  Introduction  .  1 

Background  .  1 

Objective  .  2 

Scope . 5 

Overview  .  6 

II.  The  H-function .  8 

Definition  .  8 

Convergence  Conditions  .  9 

Properties . lo 

Reciprocal . 10 

Argument  to  a  Power . lo 

Multiplication  by  an  Argument  to  a  Power.  10 

Reduction . 11 

Special  Cases  .  13 

Integral  Transforms  .  13 

Mathematical  Functions  .  13 

Statistical  Distributions  .  16 

Mellin  Transformation  and  Moments  .  23 

III.  Methods  of  Density  Estimation  .  27 

Generalized  Families  .  27 

Parametric  Estimation  Techniques  .  30 

Method  of  Moments  . . 32 

Method  of  Maximum  Likelihood  .  33 

Method  of  Least  Squares  .  35 

Other  Methods . 36 

Technique  Selection . 39 

Method  of  Maximum  Likelihood  .  39 

Method  of  Moments . 44 


Pa^e 


IV.  Nonlinear  Solution  Technique  .  45 

Fixed-Point  Iteration  .  45 

Newton's  Method  .  47 

Modified  Newton's  Method  .  50 

Steepest  Descent  Method  .  51 

Levenberg-Marquardt  Method  .  52 

Broyden's  Procedure  .  53 

Powell's  Method  .  54 

V.  Methodology  . . 56 

Generation  of  Equations  .  56 

Progr- n  Development  .  .  60 

s  ase  One  . . 60 

Phase  Two  . . 61 

Phas?>  Three  . . 66 

Phase  Four  . . 72 

Program  Limitations  .  75 

Data  Generation  .  . . 78 

Procedure . 80 

VI.  Results . 83 

Measure  of  Merit . 83 

Graph  Description . 83 

Demonstration  Runs . 89 

VII.  Conclusions  and  Remarks . 99 

Suxmary . 99 

New  Findings  . . 101 

Recommendations  for  Further  Research  .  101 


Appendix  A*  Outlines  of  Proofs  of  New  Findings  ....  103 


Appendix  B:  Advanced  Mathematical  Functions  and 
Statistical  Distributions  Expressed 
as  H-functions  . 112 

Appendix  Cs  Computer  Program  .  115 

Appendix  Dt  Comparison  of  Estimated  H-function  and 

Actual  Data . 129 

Bibliography  .  167 

Vita . 171 


iv 


Figure  Page 

1.  Program  Flowchart . 71 

2.  Two  Densities  with  Approximately  the  Same 

First  Four  Moments . 78 

3.  Sample  Graph  . . 84 

4A.  Gamma  (0=3,  0=2)  n=20 . 130 

4B.  Gamma  (9=3,  0=2)  n=150 . 131 

5A.  Exponential  (0=1/2)  n=20 . 132 

5B.  Exponential  (0=1/2)  n=150  133 

6A.  Chi-Square  (9'=4)  n=20 . 134 

6B.  Chi-Square  (0'=4)  n=150  135 

7A.  Weibull  (9=3/2,  0=1)  n=20 . 136 

7B.  Weibull  (0=3/2,  0=1)  n=150 . 137 

8A.  Rayleigh  (0=4)  n=20 . 138 

8B.  Rayleigh  (0=4)  n=150  .  139 

9A.  Maxwell  (0=2)  n=20 . 140 

9B.  Maxwell  (0=2)  n=150  141 

IOA.  Half-Normal  (9=1)  n=20 . 142 

IOB.  Half-Normal  (0=1)  n=150  143 

IIA.  Beta  (0=2,  0-3)  n=20 . 144 

IIB.  Beta  (0=2,  0=3)  n=100 . 145 

12A.  Power  Function  (9=3)  n=20 . 146 

12B.  Power  Function  (0=3)  n=100  .  147 

13A.  Uniform  n=20 . 148 

13B.  Uniform  n=100 . 149 


v 


14A.  Half-Student  (©*16)  n=20 . 150 

14B.  Half-Student  (©*16)  n*150  151 

15.  Bessel  (©*1,  0-1)  n-150 . 152 

16A.  2z2e  z  n*20 . 153 

16B.  2z2e"z  n*150 . 154 

17A  ze~z  n*20 . 155 

17B.  ze“z  n*150 . 156 

18A.  z2  n*20 . 157 

18B.  z2  n*100 . 158 

19A.  z  n-20 . 159 

19B.  z  n-100 . 160 

20A.  I  n*20 . 161 

z 

20B.  i  n*100 . 162 

z 

21A.  fz(l-z)  n-20 . 163 

21B.  /z(l-z)  n-lOO . 164 

22A.  Vz(l-z)2  n*20 . 165 

22B.  fz(l-z)2  n-lOO . 166 


vx 


rVt^  in v  viiVir  r 


List  of  Tables 


Page 


Statistical  Distributions  .  21 

Distributions  with  Constant 

H-£unction  Parameters  .  65 

Generalized  Statistical  Distributions  .  81 

Estimated  Mean  Squared  Errors  .  88 


Abstract 


, 

-"The  H-function  is  the  most  general  special  function, 
encompassing  as  specific  cases  many  mathematical  functions 
and  nearly  every  continuous  statistical  distribution  defined 
over  positive  x.  A  general  procedure  is  developed  to  esti¬ 
mate  the  parameters  of  the  H-function  which  gives  the  best 
fit  to  a  set  of  data.  The  technique  creates  a  system  of 
nonlinear  equations  from  the  method  of  moments  and  uses 
Powell's  quasi-Newton  hybrid  algorithm  to  solve  the  equa¬ 
tions.  A  computer  program,  which  can  accept  both  raw  data 
or  previously  calculated  moments,  implements  the  general 
process.  Several  new  theoretical  results  are  also 
presented . 


AN  APPLICATION  OF  THE  H-FUNCTION 
TO  CURVE-FITTING  AND  DENSITY  ESTIMATION 


I  introduction 

Analysts  cannot  effectively  work  with  raw  data.  Before 
statistical  analysis  can  be  accomplished,  data  must  be  sum¬ 
marized  in  a  convenient  form.  The  most  common  way  of 
achieving  this  is  to  fit  data  with  the  best  function  or 
statistical  distribution.  For  example,  if  the  analyst  wants 
to  model  the  time  between  arrivals  to  a  queue  for  a  computer 
simulation,  he  infers  the  true  probability  distribution  from 
a  random  sample  of  observations.  This  process  of  statis¬ 
tical  inference  usually  involves  estimating  both  the  form 
and  the  parameters  of  the  distribution.  Only  after  the  best 
form  has  been  found  can  the  analyst  begin  to  use  standard 
statistical  tools  such  as  confidence  intervals  and 
hypothesis  tests. 

Background 

Currently,  analysts  use  the  following  procedure  when 
attempting  to  fit  mathematical  functions  or  statistical 
distributions  to  data: 

1.  Plot  the  data, 

2.  Hypothesize  a  particular  type  of  function, 

3.  Estimate  the  parameters,  and 
Test  for  goodness-of-fit  using  an  appropriate  test. 
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This  procedure  has  two  major  shortcomings.  First,  the  ana¬ 
lyst  must  test  each  specific  function  separately.  For  exam¬ 
ple,  one  must  perform  two  individual  tests  if  the  data  is 
suspected  of  being  from  either  a  gamma  distribution  or  a 
half-normal  distribution.  Second,  several  distributional 
forms  and  many  combinations  of  parameters  may  fj*-  the  data. 
In  the  same  example,  if  both  goodness-of- fit  t  its  fail  to 
reject  the  hypothesized  distribution,  then  ne  her  can  be 
eliminated  as  the  true  population  distributio  Further¬ 
more,  there  is  no  chance  of  finding  the  distribution  which 
best  fits  the  data  unless  all  possible  distributions  are 
tested. 

To  solve  these  problems,  analysts  can  use  a  general 
special  function,  called  the  H-function,  because  it  includes 
many  mathematical  functions  and  statistical  distributions  as 
special  cases.  An  analyst  could  simultaneously  consider  the 
special  cases  by  simply  fitting  the  H-function  to  the  data. 
Although  this  idea  seems  logical,  the  H-function  has  never 
been  applied  to  the  current  procedure  for  curve- fitting 
because  of  the  newness  and  difficulty  of  the  H-function 
theory . 

Objective 

The  theory  has  progressed  enough  for  an  application  of 
the  H-function.  Special  theorems  allow  the  application  to 
curve-fitting  without  a  full  understanding  of  the  complex  H- 
function  theory.  Therefore,  the  overall  objective  of  this 
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research  is  to  develop  an  efficient  and  effective  method  to 
apply  the  H-function  to  the  current  procedure  for  curve¬ 
fitting  and  density  estimation.  To  reach  this  objective, 
two  main  goals  need  to  be  accomplished. 

First,  a  method  of  estimating  the  parameters  of  the  H- 
function  has  to  be  found.  The  parameters  uniquely  define 
the  H-function  and  therefore,  with  knowledge  of  the  param¬ 
eters,  the  H-function  can  be  explicitly  evaluated  and 
graphed . 

At  this  time,  a  subtlety  needs  to  be  discussed.  The 
analyst  may  also  know  that  the  data  comes  from  a  statistical 
distribution.  Because  the  method  estimates  parameters  based 
on  a  finite  number  of  data  points,  the  parameters  that 
appear  in  the  H-function  may  not  exactly  define  a  statisti¬ 
cal  distribution.  For  example,  a  parameter  required  to  be  1 
for  a  chi-square  distribution  may  be  1.01  when  estimated. 
The  other  estimates  may  also  be  slightly  off. 

However,  the  estimates  do  not  need  to  be  reevaluated. 
The  H-function  found  will  fit  the  data  better  than  any 
forced  chi-square  distribution.  If  necessary,  the  analyst 
should  think  of  the  unnamed  H-function  as  a  nearly  chi- 
square  distribution  and  proceed  with  further  analysis  given 
this  new  information. 


The  second  main  goal  is  to  determine  the  efficiency  and 
effectiveness  of  the  H-function  curve-fitting  procedure. 
This  part  of  the  research  is  somewhat  subjective  as  can  be 
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seen  from  the  following  definitions  of  efficiency  and 
effectiveness . 

Efficiency  is  measured  in  terms  of  not  only  the  number 
of  separate  tests  required,  but  also  the  difficulty  of  each 
test.  Because  the  H-function  is  a  general  function,  the 
application  of  the  H-function  reduces  the  number  of  separate 
tests  required.  Since  the  H-function  method  can  be  imple¬ 
mented  on  the  computer,  curve-fitting  will  not  be  difficult. 
Therefore,  the  new  technique  should  be  efficient. 

Effectiveness  is  measured  in  terms  of  the  ability  of 
the  method  to  find  the  “best"  function.  Analysts  have  many 
techniques  they  use  to  compare  functions  in  order  to  deter¬ 
mine  which  one  "best"  fits  the  data.  All  involve  some 
measure  of  the  error  between  the  proposed  function  and  the 
data.  Examples  include  the  absolute  distance,  the  maximum 
absolute  distance,  and  the  square  of  the  distance.  Because 
of  its  common  use,  the  estimated  mean  squared  error  will  be 
our  criterion  for  measuring  the  closeness  of  the  H-function 
to  the  data.  A  more  formal  definition  follows  in  Chapter  6. 

Three  of  the  problems  which  influence  effectiveness  are 
the  sensitivity  to  the  number  of  data  points,  the  inaccuracy 
of  higher  degree  moments,  and  the  necessity  for  an  initial 
guess  of  the  parameters.  The  size  of  the  sample  is  a  prob¬ 
lem  for  all  statistical  methods.  For  the  more  common  cases, 
only  four  to  six  moments  have  to  be  found.  The  third  prob¬ 
lem,  involving  the  initial  guess,  can  be  controlled  by 


checking  it  against  a  set  of  H-function  convergence  condi¬ 
tions.  Effectivenes.  should  increase  because  the  nonlinear 
solution  method  results  in  accurate  convergence  most  of  the 
time . 

Scope 

The  H-function  is  only  applicable  for  continuous  func¬ 
tions  defined  over  positive  values  of  x.  This  is  not  as 
serious  a  restriction  as  it  first  appears.  Methods  exist 
which  can  fold  a  symmetric  distribution  and  move  its  axis  of 
symmetry.  Therefore,  distributions  like  a  normal  or  Stu¬ 
dent's  t  can  be  evaluated  with  ‘he  H-function  as  half-normal 
or  half-student.  Such  transformations  are  not  the  subject 
of  this  thesis.  Distributions  such  as  the  half-normal  will 
only  be  analyzed  directly. 

The  H-function  also  can  be  designated  with  a  certain 
order.  Order  is  defined  as  the  sum  of  the  number  of  gamma 
terms  in  the  definition  of  the  H-function.  This  definition 
will  be  seen  later  in  Chapter  2.  For  programming  purposes, 
the  highest  order  covered  in  this  effort  is  five.  Again, 
this  is  not  a  serious  restriction.  In  fact,  most  known 
statistical  distributions  can  be  described  by  an  H-function 
with  order  one  or  two.  For  more  advanced  mathematical 


functions  such  as  arcsin  or  arctanh,  the  H-function  still 
only  needs  to  be  of  order  four. 
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Overview 


Chapter  2  contains  a  general  discussion  of  the  H- 
function  including  the  definition,  some  special  properties, 
and  many  special  cases.  A  deliberate  attempt  is  made  to 
avoid  discussion  of  as  much  of  the  complex  theory  as 
possible . 

Chapter  3  presents  the  results  of  the  curve-fitting 
literature  review.  It  concludes  with  the  selection  of  the 
appropriate  parametric  curve- fitting  method  to  use  with  the 
H-function,  the  method  of  moments. 

The  method  of  moments  produces  a  system  of  nonlinear 
equations  that  needs  to  be  solved.  The  system  is  nonlinear 
because  each  equation  involves  products  and  quotients  of 
gamma  functions  where  the  unknowns  are  in  the  arguments.  In 
Chapter  4,  a  historical  survey  of  nonlinear  solution  tech¬ 
niques  is  given.  The  technique  known  as  M.J.D.  Powell's 
hybrid  algorithm  is  selected.  This  algorithm  is  already 
available  on  the  AFIT  CDC  Cyber  750  computer  in  an  IMSL 
routine  called  ZSPOW. 

Chapter  5  discusses  the  development  of  a  computer  pro¬ 
gram  that  accomplishes  the  main  objective  of  the  research. 
The  program  estimates  the  parameters  of  the  H-function  from 
a  set  of  univariate  data  or  paired  data  after  the  data  has 
been  converted  to  moments.  Once  the  parameters  are  esti¬ 
mated,  the  H-function  can  be  explicitly  evaluated  and 
graphed.  The  H-function  is  evaluated  by  placing  the 
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estimated  parameters  into  a  program  which  has  the  capability 
of  determining  the  value  of  the  H-function  for  any  x  >  o 
[7:Appendix  B].  The  H-function  will  be  evaluated  at  the 
same  x  values  as  the  data  points.  This  will  be  the  longest 
chapter  since  it  summarizes  the  complete  methodology  of  the 
thesis  effort. 

Chapter  6  describes  graphs  which  contain  the  estimated 
H-function,  the  actual  data  points,  and  the  measure  of  the 
fit  of  the  H-function  to  the  data.  Also,  this  chapter 
summarizes  the  measure  of  fit  for  each  graph  in  tabular 
form. 

Finally,  we  reach  conclusions  about  the  efficiency  and 


effectiveness  of  the  H-function  curve-fitting  procedure. 
Also  in  the  final  chapter,  new  findings  are  highlighted  and 
further  studies  are  recommended. 


II 


The  H-function 


Definition 

The  H-function  is  the  most  general  special  function, 
encompassing  as  special  cases  most  of  the  other  special 
functions  and  elementary  functions  of  mathematics  and  nearly 
every  continuous  statistical  distribution  defined  over 
positive  x. 

The  H-function  is  defined  by  either  of  the  two  forms: 


H(z) 


m  n 

H  C  z :  [  (  , A^) } ,  i=*l, ...  ,p? 

p  q 


{ ( b  j  ,  a  j  )  }  ,  j—  l#»»*»qD 


2^r  f 

Jr 


m 


n 


TT  r(bi+B,s)  TT  Td-ai-AiS) 

j=l  J  J  i»l _  z~: 


ds 


p  q 

C1  TT  nai+Ais>  TT  Rl-N-B.s) 
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2tri  /  P  q 

Jq2  TT  r(»i-Ais)  TT  Rl-kj+B-is) 
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zs  ds 


where  z  and  all  a^  and  bj  are  real  or  complex  numbers,  all 
A^  and  Bj  are  positive  real  numbers,  and  m,  n,  p,  and  q  are 
integers  such  that  0<m<q  and  0<n<p.  Empty  products  are 
defined  to  be  unity  (1).  The  path  of  integration,  C^,  is  a 
contour  in  the  complex  s-plane  from  w-i«  to  w+i»  ,  such 

m 

that  all  Left  Half-Plane  (LHP)  poles  of  JT  Rb-i+B-iS)  lie 
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to  the  left  of  C-^  and  all  Right  Half-Plane  (RHP)  poles  of 
n 

TTr<l-ai  “  A:s)  lie  to  the  right.  Similarly,  C2  is  a 
i=l 

contour  running  from  v-ioo  to  v+i«  ,  such  that  all  RHP  poles 
m 

of  TT  r  (b-i  ”  Bjs)  lie  to  the  right  of  C2  and  all  LHP 
j=l  J  J 

n 

poles  of  fjPd-aj  +  Ajs)  lie  to  the  left  (26:2-3;  37:195; 
i=l 

7:32) . 

We  will  use  the  first  equation  in  Eq  (2.1)  as  our 
definition,  although  the  two  definitions  are  equivalent. 
When  there  is  little  chance  of  confusion,  we  will  often 
abbreviate  the  H- function  as  either 


m  n 

H  (z) 
P  q 


or 


m  n 

H  [z:{(ai,  A^};  ((bj,  B  ■ ) }  ]  . 

P  q 


We  will  define  the  order  of  an  H-function  as  p+q.  This 
represents  the  number  of  pairs  of  parameters  {(a^,  A^)}  and 
{(bj,  B j ) }  where  each  pair  represents  a  gamma  function  in 
the  integrand  in  Eq  (2.1). 

Convergence  Conditions 

For  our  purposes,  the  H-function  defined  in  Eq  (2.1)  is 
valid  if  (26:3;  7:72): 


n 

5, 


I  Ai +  Z  Bj  - 


j-W+1 


Bj  -  f  Ai  >  0 
J  i=*n+l 


and  a  Cj^  line  can  be  placed  between  the  LHP  and  RHP  poles. 
More  stringent  conditions  are  developed  later  in  Chapter  5 
when  D=0.  When  D<0,  the  H-function  is  not  defined  because 
the  infinite  sum  which  can  represent  it  does  not  converge. 
Properties 

Although  we  will  rarely  use  the  following  properties  in 
this  thesis,  we  will  state  them  for  the  reader's  future  use 
(37:196;  26:4;  7:33-34): 

Reciprocal  Property. 


H  C-:{(ai#  A^};  {(bj,  Bj)}] 

P  q 


=  H  Czs((l-bj,  Bj)};  { ( l-a±,  A±)}] 

q  p 


Argument  to  a  Power  Property. 
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P  q 
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q  P 


where  c<0 


Multiplication  by  an  Argument  to  a  Power  Property. 
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Reduction  Property. 

If  a  pair  of  "A"  terms  and  a  pair  of  "B“  terms  in  an  H- 
function  are  identical  and  one  is  in  the  numerator  and  the 
other  is  in  the  denominator,  then  it  is  equivalent  to  an  H- 
function  with  a  lower  order.  Specifically  (26:4;  7:34-35): 


m  n 

H  C  z :  {  (  a^ ,  Aj. )  }  ;  •••#  (b__  ^,B__  ^)#  (a^,A^)3 

p  q 

m  n-1 

*  H  [z:  (a2#A2 )  #  •  ••»  (s^.A—);  {(b-i»B.i)}3 

p-1  q— 1  Z  *  P  P  3  3 


provided  n>0  and  q>m. 


m  n 

H  Cz:  (a^,A^) ,  ...»  ( #  Ap_j_ ) »  ;  {(Oj,Bj)}J 


m-1  n 

*  B  [z:  { (  ai#  Aj_) } ;  (b2»B2^»  •••#  (b_,Bg)3 

p-1  q-1 


provided  m>0  and  p>n. 

We  also  found  another  way  in  which  the  H-function  can 
reduce  to  one  of  lower  order.  If  any  A^  or  Bj  is  close 
enough  to  zero,  that  gamma  term  in  the  integrand  of  Eq  (2.1) 
is  essentially  a  constant.  Thus, 
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■  •  ■  'V'  ^ 


H  [z]  ~P(to1)  H  [z:  {(ai,Ai)}  ;  (b2,B2)f 

P  q  P  q-1 


•  •  •  *  ( bq , Bq ) ] 


for  B-^-O  and  m>l. 


m  n  m  n 

H  [z]  ~  ■  ■  ■■■  H  C z :  { (  a^ / 

P  q  I  U-Dq'  P  q-1 


•••»  (bq_i#Bq_i)3 


for  Bq~0  and  m<q. 


m  n  m  n-1 

H  [z]  ~  P*  ( 1— a^ )  H  [z:  (a2»A2)  ,  •  •»,  (a_,Ap)7 

p  q  p-i  q 


Ubj.Bj)}] 


for  A^~0  and  n>l, 


m  n 


H  CzQ  ~  -r  ■  H  [z:(a^,A^)/  •  ••#  (  3p_^  #  Ap— i )  7 

p  q  I  'apJ  p-1  q 

{ (bj ,Bj ) } ] 

for  Ap~0  and  n<p. 

Therefore,  there  are  two  ways  in  which  the  H-function 
can  reduce  to  a  lower  order.  Gamma  terms  in  the  integrand 
of  Eq  (2.1)  could  cancel  in  the  numerator  and  denominator  or 
they  could  reduce  to  constants.  These  reduction  properties 
could  be  useful  in  allowing  a  less  restrictive  assumption  of 
the  values  of  m,  n,  p,  and  q  when  fitting  the  H-function  to 


Integral  Transforms. 


The  Laplace  (Lr)  and  Fourier  (Ft)  transforms  of  an  H- 
function  are  also  H-functions  (37*199-201;  7s35): 


Lr  H  (cz) 

(  pq  ) 

m  n+1 

=  —  H  [£: (0,1), {(a,,  A.:)};  {(bi#  Bi ) } ] 

r  p+1  q  r  J 


Using  the  property  of  multiplication  by  an  argument  to  a 
power,  this  becomes: 


m  n+1 


*  i  H 


p+1  q 


[£:  (1,1),  {(aj+Ai,  A^};  {(bj+Bj,  Bj)}] 


By  the  reciprocal  property,  this  can  be  rewritten  as : 


n+1  ra 

«±H  [£:  { ( 1-b.i-B.i ,  Bi)};  (0,1),  {(l-a^-Ai,  A.)}1 

c  q  p+1  c  3  3  3  X  l  1 

Im  n  I 

H  (cz)  } 


n+1  m 

Ih  C— :((l-bi-Bi,  B,)};  (0,1),  {(l-aj-A,,  A±) }  ] 

c  q  p+1  c  3  3  3  111 


Mathematical  Functions. 

The  following  elementary  mathematical  functions  can  be 
expressed  as  H-functions  (26:10,151-152;  7:39-41,124): 


z°e~z  -  H  [zs ; (b, 1 ) ] 
O  1 


H  B  10 

L  zB  e"z  =  H  [z:;{b#B)] 


zb  =  H  Cz:(b+l,l);  (b,  1 ) ] 

1  1 

1  0 

zb  (l-z)+a  =  r(a+l)  H  [z : ( a+b+1, 1 ) ;  (b,l)] 

1  1 


zb  (l+z)“a  -  H  [z: (b-a+1,1);  (b,l)] 

1  1 


sin (z)  -  H  [Z:;(I,±),  (0,1)] 
2  0  2  1  1  1  1 


sinh(z)  »  H  (0,1)] 

2  0  2  2  2  2  2 


cos(z)  »  H  Czs;(0,I),  (i,I)] 

2  o  2  2  2  22 


arcsin(z) 


1  0 

H  [i #: 

7  (0,1),  (1,1)] 

0  2  2 

2  2  2 

12  .  . 

•  _*  H 

4*^  2  2 

[izs (l,i)  ,  (l,i); 
2  2 

(i.i: 
2  2 

12  _  _ 

-±-  H 

[zs(l,4),  (1,1); 

(±.±> 

4V?C  2  2 

2  2 

2  2 

arctan(z)  ■  I  H  [z:(l,i),  (i,i)  ;  (i,i) ,  (0,-i.)] 


1  2 

arctanh(z)  =  -i  H  [izs(l,A),  (i,i)  ;  (i,!)  ,  (0,1)] 
4  2  2  2  2  2  22  2 

1  2 

log(l±z)  =  H  [-z: (1,1),  (1,1);  (1,1),  (0,1)3 
2  2 


12  0 

-H  [z: (1,1),  (1,1);  (0,1),  (0,1)]  0<z<l 

2  2 

0  2 

H  [zs (1,1),  (1,1);  (0,1),  (0,1)]  z>l 

2  2 


It  should  be  noted  that  the  formulas  given  above  are 
not  consistent  with  those  in  Mathai  and  Saxena  for  arc- 
sin(z),  arctanh(z),  and  log(l-z).  There  were  apparently 
typographical  errors  in  their  book.  We  believe  that  the 
correct  formulas  are  given  above.  We  verified  these  as  the 
correct  formulas  by  summing  the  residues  of  the  H-functions 
given  above.  In  each  case  we  obtained  the  correct  infinite 
series.  The  outlines  of  these  proofs  are  provided  in 
Appendix  A. 

We  were  able  to  prove  a  generalization  to  the  loga¬ 
rithmic  function  log(z).  It  is  given  below  and  the  proof  is 
outlined  in  Appendix  A. 


log(z)  *  -u2  H  Czs(l,u),  ( l»u) ;  (0,u) ,  (0, u) ]  0<z<l 

2  2 


u2  H  Cz:(l,u),  ( 1, u) ;  (0,u),  (0,u)] 
2  2 
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A  similar  generalization  applies  to  the  power  function 
z*3.  it  is  stated  below  and  the  proof  is  a  special  case  of 
the  proof  for  the  generalization  of  the  Power  Function 
probability  density  function  (p.d.f.)  in  Appendix  A. 


zb  =  uH  [z:(ub+l,u);  (ub,u)3 
1  1 


The  H-function  also  includes  as  special  cases  many 
advanced  mathematical  functions.  Since  these  will  not  be 
dealt  with  in  this  thesis,  they  are  provided  for  the 
reader’s  benefit  in  Appendix  B. 

Statistical  Distributions. 

Consider  a  continuous  random  variable  X  whose 
probability  density  function  (p.d.f.)  is  given  by 


f  v  ( x )  - 


m  n 


k  H  ( cx ) 

P  q 


cx  C  S 


0  otherwise 

r00 

where  k  and  c  are  constants  such  that  /  fx(x) 


dx  =  1  and 


-00 


m  n 

S  is  a  subset  of  the  positive  real  line  for  which  H  (cx) 

P  q 


is  convergent.  Then  the  random  variable  X  is  said  to  be  an 
H-function  variate  or  a  random  variable  with  an  H-function 
distribution  (37s 200;  7:84). 

Many  common  statistical  distributions  are  special  cases 


of  the  H-function  distribution.  These  include  (37:164,  202- 

207;  7:85-87,93-94): 


Gamma  p . d . f . 


f  ( x  I  © ,  0 )  = 


_  0  x®-1  e”0x 


r<©> 


=  -1 —  H  [0x: ; (0-1, 1)] 

H©)  O  1 

Exponential  p.d.f.  (Gamma  p.d.f.  with  ©  =  1) 

(Weibull  p.d.f.  with  0*1) 

f  ( x I  0 )  =  0e”0x 
1  O 

=  0H  [0x: ; (0,1)3 
O  1 


x  >  O 


©,0  >  O 


x  >  O 


0  >  O 


Chi-Square  p.d.f.  (Gamma  p.d.f.  with  0  =  £  and  0  =  i) 


©'  ©'-1  -x 
f(xle')  *  [22  rt®')]-1  x2  e  2 


=  C 2  nt')!'1  H  Cx*;(|' -1,1)3 


2  2 


x  >  O 


0 1  >0 


Weibull  p.d.f. 


f  (  X  I  0 , 0  )  * 


©Gx9”1  e"0x 

I  1  o  1 


=  0®  H 


[09x: ; (1-1,1) 3 


x  >  O 


0,0  >  O 


Rayleigh  p.d.f.  (Weibull  p.d.f.  with  ©  =  2) 

_  2 

f(x|0)  =  20x  e~0x 
1  0 

=  V0  H  [V0xs;(4,l)] 

0  1  22 


Maxwell  p.d.f. 


f(xle) 


93VtF 


2 

©Vi? 


1  0 
H 

0  1 


cfwu.in 


x  >  0 
0  >  O 


x  >  O 

©  >  O 


Half -Normal  p.d.f. 


f  ( x | ©) 


.  X2 

2  e  2©2 
V2tt© 


1  O 

53  -~=-  H  C — — x:  ;  (0,-i.)  ] 

Vine  O  1  V2©  2 


X  >  O 

©  >  o 


Beta  p.d.f.  of  the  first  kind 


f(xle,0)  -  ^(8+0^  *8-1  u.x)0-i 


f<e) 


l  o 
i 

l  l 


=  .pQ-tll  H  [x:  (©+0-1,1);  (e-lfl)] 


0  <  x  <  1 


©,  0  >  O 


Power  Function  p.d.f.  (Beta  p.d.f.  with  0=1) 


f(x|6)  = 

Uniform  p.d.f 

f(x)  =  1 

Half-Cauchy  p 
f(x|e)  = 

Half-Student 

f(xl©)  = 


ex®"1 


1  o 

9H  [x:(9,l);  (0-1,1)]  O  <  x  <  1 

1  1 

0  >  O 


.  (Beta  p.d.f.  with  0=0=1) 

(Power  Function  p.d.f.  with  ©  =  1) 


1  O 

H  [x:(l,l)y  (0,1)]  O  <  x  <  1 

1  1 


-d.f . 
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/n-(02+x2) 

1  1 

~  H  C^:(0,i);  (0,i)]  x  >  O 

©ff-  x  0  2  2 

0  >  O 


p.d.f. 


2 


2 

VeS?  f(|)  (i  +  2S_) 

&  b 


0+1 

2 


1  1 


H 


1  1 


r_x. i i-e 

Ve  2 


•¥ 


(O,1)] 


x  >  O 
9  >  0 


F  p.d.f 


p .  0+0 .  |  %  f-1 

et  \  r.\  ^  2  ^  82  02  X2 

£(xle,0)  =  - ± - -  — — ^ 

rtf)  n§>  (ex+0)2 


- JL  — -  H  [®*:(-f,l);  (|-1#1)3 

0f(|)p(|)  1  1  0  2  2 


x  >  0 


0,  0  >  O 


Beta  p.d.f.  of  the  second  kind 


f  (x|e,0)  =  (0/0)e  — xftz1.  _ 

n©)R0)  u+|^> 


__J! -  H  [g*:(-0,l)y  (e-1,1)]'  X  >  o 

ep(8)p(0)  l  l  e 

©,  0  >  o 


These  results  are  summarized  in  Table  I  (8:300). 

We  were  able  to  prove  a  generalization  to  the  above 
formulas  for  the  Power  Function  p.d.f.  and  the  Uniform 
p.d.f.  This  generalization  was  not  previously  known  and  the 
proofs  are  outlined  in  Appendix  A. 


Power  Function  p.d.f. 


f(x|e)  = 


=  U0H  [x:(u(©-l)+l,u);  (u(©-l),u)]  0  <  x  <  1 

1  1 

0,  u  >  0 


VVAV.'V.iSSV 


VVV-.VV.ViV 


f  1111  e/[0P(e/2)p(0/2)]  e/0 

Beta  (second  kind)  1111  0/[©P(8)  P(0) ]  0/e 


Uniform  p.d.f.  (Power  Function  p.d.f.  with  8=1) 


f(x)  =  1 


1  0 

UH  Cxs(l,u);  (0 , u) ] 
1  1 


0  <  x  <  1 


u  >  0 


We  were  also  able  to  represent  the  Pareto  p.d.f.  as  an 
H-function.  This  is  another  new  finding  not  previously 
known.  Again,  the  proof  is  outlined  in  Appendix  A. 


Pareto  p.d.f. 


f(x|e)  =  ex-0"1 


=  ©H  [x: (-9, 1) ?  (-0-1,1)] 

1  1 


x  >  1 


8  >  0 


This  p.d.f.  can  also  be  generalized  as  above.  The 
proof  of  the  result  is  given  in  Appendix  A. 


Pareto  p.d.f. 


f ( x I o )  =  ex-9"1 


=  U0H  [x: ( l-u( 1+0) ,u) ;  (-u(l+0),u)]  x  >  1 
1  1 

0,  u  >  o 


The  Bessell  p.d.f.  and  General  Hypergeometric  p.d.f. 
can  also  be  expressed  as  H-functions.  These  are  listed  in 
Appendix  B  since  we  won't  use  them  in  this  thesis.  However, 
we  will  use  the  Bessel  p.d.f.  to  verify  and  validate  our 
computer  program. 


xV 


If  a  random  variable  has  a  p.d.f.  that  can  be  expressed 


as  an  H-function  distribution,  then  its  cumulative  distribu¬ 


tion  function  (c.d.f.)  can  be  easily  written  as  one  minus 


another  H-function.  Specifically,  if 


f(x)  =  k  H  [cx:{(ai#  Ai)};  {(bj,  B j  ) }  , 

p  q 


then  (37:243;  7:102-107): 


F  ( X )  as  1  —  is  H 


m+1  n 


p+1  q+1 


[cx^Uj+Ai,  Ai)),  (1,1);  (0,1) 


{(bj+Bj,  Bj)}] 


Another  powerful  characteristic  of  H-function  variates 


is  that  products,  quotients,  and  rational  powers  of  indepen¬ 


dent  H-function  variates  are  also  random  variables  with  H- 


function  distributions.  The  exact  formulas  are  lengthy  and 


will  not  be  presented  here.  The  reader  is  referred  to 


Springer  (37:207-219)  or  Cook  (7:90-92).  In  his  disserta¬ 


tion,  Cook  presented  a  computer  program  that  will  evaluate 


these  combinations  and  graph  the  p.d.f.  and  c.d.f.  of  the 


resulting  H-function  distribution. 


Mellin  Transformation  and  Moments 


The  Mellin  transform  (M„)  of  an  H-function  is  given  as 


(37:198-199;  7:35): 
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*11  11 

TT  ftbi  +  B.a)  TT  TU-ai  -  A^s) 

Mg{H(cx) }  »  j-1  3  3  i-1 _ 1  1  c-s 

p  _  q 

TTnai  +  Ajs)  tt  Hi-b,  -  BjS) 

i=n+l  j=m+l 


For  continuous  random  variables  defined  over  positive  x,  the 
moments  about  the  origin  are  (37:201-202;  7:29,108-109): 


/lr  *  E(Xr)  =  /  xr f ( x ) i 


Mr+i^ f(x) ) 


Therefore,  the  moments  of  the  H-function  distribution  are: 


ftj.  *  Mr+1  {kH  (cx)}  =  kMr+1  {H  (cx)}  - 


-  Ji— I(r+1) 


r— 0 1 1#  •  •  • 


where 


AM  A* 

TT  Rbj  +  Bjs>  TT  H1"3!  -  Ais> 

I(s)  -  j-1  J  3  i-1 _ 

p  q 

TT  n«i.  +  Ais>  TT  Rl-bj  -  B^s) 
i=n+l  j=m+l  J  J 

The  above  formula  is  especially  useful  for  finding  the 
constant  k  in  an  H-function  distribution.  Using  the  fact 
that  the  zeroth  moment  of  a  statistical  distribution  must 
equal  one, 

m  n 

JJ„  *  1  *  M,  {kH  (cx)}  *  *1(1) 


or  (7:109): 


p  _  q 

TT  Rai  +  A±)  JT  R !-»> j  -  Bj) 

c  «  c  i*n+l _  j=m+l _ _ _ 

I  ( 1 )  m  n 

tt  rib,  +  b,)  tt  ru“ai  -  Aj_) 

j=l  i=l 


The  moments  of  H-function  distributions  carry  their 
usual  statistical  meaning.  For  example,  the  zeroth  moment 
is  the  area  under  the  p.d.f.  over  the  appropriate  range. 
The  first  moment  is  the  mean  of  the  distribution.  The 
second  moment  is  a  function  of  the  mean  and  variance  of  the 
distribution. 

We  could  also  define  moments  of  a  mathematical  function 
f(x)  defined  over  positive  x  as: 

Alr  =  /  xr  f ( x )  dx  (2.2) 

o 

When  the  function  can  be  represented  as  an  H- function,  these 
moments  could  be  found  by  the  same  formula  as  for  statisti¬ 
cal  distributions.  However,  they  would  not  always  have  the 
statistical  meaning.  The  first  moment  of  a  function  could 
not  always  be  interpreted  as  the  mean  or  balance  point  of 
that  function.  The  definition  in  Eq  (2.2)  will  be  useful  in 
Chapter  5  when  we  fit  mathematical  functions  to  data  using 
the  method  of  moments. 

This  concludes  our  discussion  of  the  definition,  prop¬ 
erties,  special  cases,  and  moments  of  the  H-function.  We 


next  examine  procedures  to  estimate  the  parameters  of  the  H- 
function.  Because  the  H-function  includes  both  mathematical 
functions  and  statistical  distributions  as  special  cases,  we 
needed  to  find  a  method  of  estimating  the  parameters  of  the 
H- function  that  would  be  applicable  for  both  curve- fitting 
and  density  estimation.  Nearly  every  method  we  found  was 
formulated  for  density  estimation,  although  some  could  also 
be  applied  to  curve-fitting. 

It  appears  that  when  trying  to  fit  mathematical  func¬ 
tions  to  data,  analysts  usually  plot  the  data  points  and 
hope  to  find  a  pattern  in  the  points  that  is  recognizable  as 
a  special  mathematical  function.  If  the  points  do  not 
exhibit  a  pattern,  the  analysts  could  still  approximate  the 
n  data  points  with  a  polynomial  of  order  0  to  order  n-1. 

Because  of  the  lack  of  techniques  for  fitting  mathema¬ 
tical  functions  to  data  and  the  abundance  of  techniques  for 
density  estimation,  we  concentrated  our  literature  review  on 
methods  of  density  estimation.  As  noted  previously,  in  some 
cases,  these  methods  may  also  be  used  to  fit  mathematical 


functions  to  data. 


Methods  of  Density  Estimation 


Methods  of  probability  density  estimation  can  be  gener¬ 
ally  classified  as  parametric  or  nonparametric.  Parametric 
density  estimation  techniques  usually  assume  that  the  form 
of  the  distribution  from  which  the  data  were  taken  is  known. 
Alternatively,  they  may  estimate  the  form  or  class  of  the 
probability  density  function  (p.d.f.)  from  the  data.  In 
contrast,  nonparametric  density  estimation  techniques  are 
not  concerned  with  the  form  of  the  distribution,  before  or 
after  the  data  are  taken. 

Generalized  Families 

Although  parametric  approaches  generally  require  an 
assumption  of  the  form  of  the  unknown  p.d.f.,  that  assump¬ 
tion  is  not  always  as  restrictive  as  it  appears.  There  are 
several  generalized  families  which  include  many  density 
functions  as  special  cases.  Assuming  that  the  sample  came 
from  a  generalized  family  is  not  nearly  as  restrictive  as 
assuming  it  came  from  a  particular  distribution. 

The  emphasis  in  the  generalized  approach  is  to  leave 
the  functional  form  of  the  unknown  density  as  unspecified  as 
possible  and  allow  the  data  to  indicate  which  special  case 
of  the  family  gives  the  best  fit  (41:13).  This  approach 
permits  model  and  parameter  estimation  to  be  considered 
simultaneously.  This  seems  to  be  beneficial  because  of  the 
close  relationship  between  a  model  and  its  parameters 


One  of  the  best-known  generalized  families  was  proposed 
by  Karl  Pearson.  This  family,  as  it  turned  out,  included 
many  of  the  more  common  continuous  univariate  probability 


densities  as  members  (41:5). 

Pearson's  system  of  frequency  curves  is  generated  by 
solutions  to  the  differential  equation 

y '  =  X+a  y 

t>0+blx+b2x^ 

where  a  and  the  b's  are  constants  (9:248-249?  37:255).  The 
system  consists  of  12  types  of  curves  and  a  set  of  rules  for 
determining  which  curve  best  fits  the  data  based  on  the 
first  four  moments  (37:255).  A  fairly  detailed  development 
of  the  system  is  given  by  Elderton  (13:38-127). 

Special  cases  of  the  Pearson  system  of  curves  include 
the  normal  distribution,  the  chi-square  distribution.  Stu¬ 
dent's  t  distribution,  the  beta  distributions  (first  and 
second  kinds),  and  the  Pareto  distribution  (9:249;  41:8). 
The  gamma  distribution  can  also  be  obtained  after  shifting 
the  origin  and  making  a  transformation  (41:9). 

Another  generalized  family  is  the  generalized  gamma 
family  defined  by  Stacy  (38).  It  consists  of  probability 
density  functions  of  the  form: 

J2_ 

(a?  xd-l  e-(x/a)P 


f  ( x | a ,d, p)  = 


where  a,  d,  and  p  are  parameters.  This  family  includes  the 
gamma,  Weibull,  Maxwell,  and  standard  normal  distributions 


as  special  cases  (38:1187).  Of  course,  the  exponential  and 
chi-squared  distributions  are  also  included  in  the  family 
since  they  are  special  cases  of  the  gamma  distribution. 
Similarly,  the  Rayleigh  distribution,  a  special  case  of  the 
Weibull  distribution,  is  also  in  the  generalized  gamma  fam¬ 
ily.  This  family  is  equivalent  to  the  first  order  H- 

1  0 

function  distribution  kH  (cx). 

0  1 

Recently,  Ramberg  et  al  (33)  proposed  a  four-parameter 
probability  distribution  whose  percentile  function  R(p)  is 
based  on  the  generalization  of  Tukey's  lambda  function: 

X3  x4 

R(p)  =  X.  +  - "  DrJg) - 1  0<  p  <  1 

x2 

The  density  functions  related  to  this  percentile  function 
can  take  on  a  variety  of  shapes,  depending  on  the  values  of 
the  X^.  This  distribution  can  represent,  or  at  least 
approximate,  the  gamma,  Weibull,  normal,  log-normal,  and 
Student's  t  distributions  (33:203,  206).  The  proposed  dis¬ 

tribution  yields  a  good  approximation  to  the  data  using  the 
first  four  moments. 

Clearly,  the  H-function  is  also  a  generalized  family 
since  it  includes  as  special  cases  nearly  every  named  con¬ 
tinuous  probability  density  defined  over  positive  x.  The 


only  named  continuous  densities  that  have  not  been  shown  to 
be  H-functions  are  the  log-normal  and  logistic  distribu¬ 
tions.  On  the  other  hand,  no  one  has  been  able  to  show  that 
these  distributions  are  not  H-functions. 

As  seen  in  Chapter  2  and  Appendix  B,  the  H-function 
also  includes  many  named  functions  of  mathematics  as  special 
cases.  Thus,  the  assumption  that  the  data  came  from  an  H- 
function  could  hardly  be  considered  restrictive. 

Still,  that  assumption  makes  fitting  an  H-function  to 
data  a  parametric  procedure.  In  this  case,  the  assumed  form 
of  the  distribution  is  an  H-function  with  particular  values 
for  m,  n,  p,  and  q.  Estimating  the  parameters  of  that  H- 
function  would  simultaneously  consider  all  of  the  densities 
which  are  special  cases  of  the  H-function  and  determine  the 
H-function  which  best  fits  the  data. 

Parametric  Estimation  Techniques 

Since  the  form  of  the  distribution  is  assumed  to  be 
known  in  parametric  density  estimation  approaches,  the  prob¬ 
lem  reduces  to  finding  point  estimates  for  the  parameters  of 
the  p.d.f.  Before  discussing  the  various  approaches  to 
parameter  estimation,  however,  we  need  to  define  certain 
desirable  properties  of  estimators. 

A 

An  estimate,  ©,  of  an  unknown  parameter,  9,  is  said  to 

A 

be  unbiased  if  E(0)=©  for  all  9.  "This  implies  that  the 


sampling  distribution  of  ©  is  centered  at  the  parameter  ©. 


That  is,  an  unbiased  estimator  ©  is  equal  to  ©  on  the 
average"  (3:388). 

A 

If  the  estimate  ©  converges  in  probability  to  ©,  it  is 

said  to  be  consistent.  Formally,  if  lim  Pr { | ©-© | < £ } =1  for 

n-»oo 

A 

any  £>0,  then  0  is  a  consistent  estimate  of  ©  (18:235-236). 

Consistent  estimators  are  not  necessarily  unbiased  and 
unbiased  estimators  are  not  necessarily  consistent.  Thus, 
neither  property  implies  the  other.  But  a  consistent  esti¬ 
mator  with  a  finite  mean  value  must  tend  to  be  unbiased  in 
large  samples  (23:5). 

An  efficient  estimator  of  ©  is  an  unbiased  estimator 
with  minimum  variance  among  all  unbiased  estimators.  A 

A 

measure  of  the  efficiency  of  the  estimator  ©a  is 

A 

efficiency  =  ^ 

v(ea) 

A 

where  V(e)  is  the  minimum  variance  of  all  unbiased  esti- 

A  A 

mators  and  v(0a)  is  the  variance  of  ©a»  "An  efficient 
estimator  is  sometimes  called  a  minimum  variance  unbiased 
estimator"  (3:388).  An  asymptotically  efficient  estimator 
is  an  estimator  that  becomes  efficient  as  the  sample  size 
increases  to  infinity. 

With  these  properties  defined,  we  can  proceed  to 


discuss  the  various  approaches  to  parameter  estimation. 


/  f 


A 


Method  of  Moments. 

The  method  of  moments  was  proposed  by  Karl  Pearson  to 
approximate  data  with  a  curve.  This  method  involves 
equating  the  moments  of  the  data  with  the  moments  of  the 
distribution,  creating  as  many  equations  as  there  are  param¬ 
eters  to  be  estimated.  The  estimates  are  then  obtained  by 
solving  these  equations  for  the  parameters  (9:497-498; 
13:12-37;  29:274-276).  "This  method  often  leads  to 

comparatively  simple  computations  in  practice"  (9:497). 

However,  one  limitation  of  the  method  of  moments  is  the 
unstable  nature  of  the  higher  moments  calculated  from  the 
data.  If  there  are  many  unknown  parameters,  then  higher 


order  moments  will  be  required  to  solve  for  the  parameters. 


Karl  Pearson  has  shown  that  "we  might 
easily  on  a  random  sample  reach  a  7th  or 
8th  moment  having  half  or  double  the 
value  it  actually  has  in  the  general 
population.  Constants  based  on  these 
high  moments  will  be  practically  idl'j. 
They  may  enable  us  to  describe  closely 
an  individual  random  sample  but  no  safe 
argument  can  be  drawn  from  this 
individual  sample  as  to  the  general 
population  at  large,  at  any  rate  so  far 
as  the  argument  is  based  on  the 
constants  depending  on  these  high 
moments"  (13:44). 


This  limitation  led  to  the  development  of  equations  to 
correct  the  raw  or  grouped  moments  (9:360-362). 


Another  concern  with  the  method  of  moments  is  the 
question  of  whether  a  finite  number  of  moments  can  uniquely 
determine  the  distribution.  Although  there  is  a  one-to-one 


correspondence  between  the  moment  generating  function  and 
the  distribution,  unless  the  moment  generating  function  is 
known,  the  moments,  in  general,  do  not  uniquely  determine 
the  distribution  function.  This  concern  is  referred  to  as 
the  problem  of  moments  (29:81). 

Although  estimates  obtained  with  the  method  of  moments 
are  sometimes  biased,  we  can  often  remove  the  bias  with  a 
simple  correction  and  thus  obtain  an  unbiased  estimate 
(9:498).  "In  general,  these  estimates  are  consistent" 
(3:389).  However,  the  asymptotic  efficiency  of  the  esti¬ 
mates  is  often  considerably  less  than  1,  which  implies  that 
they  are  not  the  "best"  possible  estimates  from  the 
efficiency  point  of  view  (9:498). 

Method  of  Maximum  Likelihood. 

The  concept  of  maximum  likelihood  was  first  introduced 
by  R.  A.  Fisher  in  1912  and  applied  to  parameter  estimation 
in  1921  (9:498). 

There  are  two  crucial  assumptions  of  the  method  of 
maximum  likelihood.  First,  the  correct  form  of  the  equation 
must  be  known  or  assumed.  Second,  the  data  must  be  a  repre¬ 
sentative  sample  from  the  whole  range  of  situations  about 
which  the  analyst  wishes  to  generalize  (10:7).  The  second 
assumption  is  common  to  all  methods  of  density  estimation. 
The  first  assumption,  however,  is  characteristic  only  of 
parametric  estimation  techniques,  although  some  parametric 
techniques  allow  a  less  restrictive  assumption  of  the  form 


of  the  equation  than  others.  This  point  was  emphasized 
earlier  in  this  chapter. 

The  method  of  maximum  likelihood  consists  of  deter- 

A  A  A 

mining  the  values  ©^,9 2'"*'ek  maximize  the  likelihood 

function  with  respect  to  ©^,...,9^.  The  likelihood  function 

is  defined  as: 

* 


n 

=  TT 

i=l 


f  (xi 


where  f  ( x^  1 0^ ,  • 9^)  is  the  p.d.f.  of  X^.  This  is  equiv¬ 
alent  to  maximizing  log  L  with  respect  to  9^...,©^  since 
both  L  and  log  L  are  maximized  at  the  same  value.  This  is 
useful  since  log  L  is  sometimes  easier  to  maximize  than  L. 

To  maximize  the  likelihood  function,  an  analyst  usually 
differentiates  L  (or  log  L)  with  respect  to  each  of  the 
unknown  parameters  ©^,...,9^.  These  derivatives  are  then 
set  equal  to  zero  and  the  resulting  system  of  equations  is 
solved  for  9^...,©^.  The  solutions  to  these  equations  are 
the  maximum  likelihood  estimators  (23:35-74;  27:183-186). 
The  primary  difficulty  with  this  method  is  that  the  system 
of  equations  often  cannot  be  solved  directly  and  the 
constants  have  to  be  found  by  numerical  approximation 
(13:252)  . 

Still,  under  some  general  conditions,  maximum  likeli¬ 
hood  estimators  are  consistent,  asymptotically  normal,  and 
asymptotically  efficient  (3:389).  Although  the  estimates 


34 


are  not  necessarily  unbiased,  many  times  they  can  be  modi¬ 


fied  so  that  they  become  unbiased  (28:186).  For  these 
reasons,  the  method  of  maximum  likelihood  is  the  most  widely 
used  density  estimation  technique  (41:13). 

Method  of  Least  Squares . 

Another  popular  technique  for  fitting  curves  to  data  is 
the  method  of  least  squares.  This  method  involves  finding 
the  constants  of  the  assumed  equation  which  minimize  the 
square  of  the  differences  between  the  actual  data  values  and 
the  values  predicted  by  the  equation  (23:75-91;  29:482-502). 

Linear  least  squares  is  a  well-developed  technique  that 
can  be  applied  to  any  form  of  equation  that  can  be  refor¬ 
mulated  through  transformations  into  another  equation  that 
is  linear  in  its  coefficients.  Daniel  and  Wood  (10:19-23) 
suggest  appropriate  transformations  to  transform  quite  a  few 
nonlinear  equations  into  a  model  that  is  linear  in  its 
coefficients . 

Nonlinear  least  squares  estimation  is  a  relatively  new 
area  developed  to  accomodate  models  which  cannot  be  made 
linear  through  transformations.  Several  methods  are  avail¬ 
able  which  use  numerical  techniques  such  as  Gauss-Newton  or 
steepest-descent  to  converge  on  a  solution  (10:9-10). 

For  the  simple  linear  model  and  when  general  assump¬ 
tions  are  made,  the  Gauss-Markov  theorem  states  that  the 
least  squares  estimators  are  the  best  (i.e.  minimum  vari¬ 
ance)  linear  unbiased  estimators  of  the  unknown  coefficients 
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in  the  model.  Further,  when  the  random  errors  in  the  model 
are  normally  distributed,  the  least  squares  estimates  are 
maximum  likelihood  estimates  and  are  of  minimum  possible 
variance  (10:7). 

Other  Methods . 

In  the  minimum  chi-square  method,  the  observations  are 
grouped  into  c  intervals  and  the  values  ©-^,...,0^  are  found 
which  minimize 

c 

■y  2  =  £  [ni  -  npi(©1, . . . ,ek) ]2 

1=1  rip.  (©  ,  . . .  ,0,  ) 
i  1  k 

where  n^  is  the  actual  number  of  observations  in  interval  i 
and  np^(01,  ...,9k)  is  the  predicted  number  of  observations 
in  the  interval,  regarded  as  a  function  of  ©lf...,©k 
(3:389).  The  asymptotic  properties  of  minimum  chi-square 
estimators  are  similar  to  those  of  maximum  likelihood  esti¬ 
mators  (23:93).  But  as  with  some  previously  mentioned  tech¬ 
niques,  the  equations  are  usually  too  difficult  to  be  solved 
analytically  and  a  numerical  technique  must  be  used.  Fur¬ 
ther,  the  observations  must  be  grouped,  even  when  dealing 
with  a  continuous  distribution,  and  it  seems  rather  wasteful 
to  impose  an  otherwise  unnecessary  grouping  for  estimation 
purposes  (23:93). 

In  Bayesian  statistics,  0  is  not  regarded  as  an  unknown 


constant,  but  as  a  random  variable. 


Thus,  it  has  a 


probability  density  function,  although  this  p.d.f.  is 
unknown.  The  objective  in  Bayesian  estimation  is  to  combine 
any  prior  information  about  the  distribution  of  e  with  the 


random  sample  before  estimating  ©.  Bayes’  method  is  well 
formulated  for  a  single  ©.  It  involves  multivariate  dis¬ 
tributions  of  ©^,...,©^  when  a  vector  of  parameters  is 
considered . 


If  g(©)  is  the  p.d.f.  of  the  parameter  which  expresses 
the  prior  information  about  ©  and 


f(x1# . . . ,xnle)  =  jj  f(xil©) 

i=l 


is  the  joint  p.d.f.  of  X 


1  i  •  •  • 


Xn,  given  ©,  then  the 


posterior  p.d.f.  of  ©,  given  the  random  sample,  is 


h(©| x, , . . . , x  )  = 
1  n 


g(©)  f . . . ,xnl©) 

'g  ( e )  f ( xx , . . . , xn | © )  d© 


where  the  integration  is  performed  over  the  possible  values 
of  ©.  The  posterior  p.d.f.  of  ©  represents  the  current 
knowledge  about  9,  incorporating  the  prior  p.d.f.  of  ©  and 
the  random  sample.  Any  measure  of  centrality  of  the  poste¬ 
rior  p.d.f.,  such  as  the  mean,  median,  or  mode,  can  be  used 
as  a  point  estimate  of  ©  (29:339-351). 

Another  relatively  simple  way  to  estimate  the  form  of 
the  distribution  is  to  use  the  graphical  method.  With  this 
approach,  the  points  of  the  empirical  (i.e.  sample)  cumu¬ 
lative  distribution  function  (c.d.f.)  are  plotted  on 


probability  paper  of  the  assumed  type  of  distribution.  If 
the  points  lie  roughly  in  a  straight  line,  then  the  correct 
form  of  distribution  was  assumed. 

With  some  types  of  probability  paper,  a  scale  is  pro¬ 
vided  to  estimate  the  parameters  of  the  p.d.f.  (22:295-308). 
Alternatively,  parameter  estimates  can  be  obtained  using 
other  techniques  such  as  maximum  likelihood. 

A  lesser-known  technique  of  density  estimation  is  the 
minimum-distance  method.  Given  a  distance  function  d(F,G) 
which  measures  how  "far  apart"  two  cumulative  distribution 
functions  F  and  G  are,  the  minimum-distance  estimate  of  9  is 
the  value  of  ©  which  minimizes  d(F(x  |e),Fn(x)),  where  Fn(x) 
is  the  empirical  c.d.f.  Although  intuitively  appealing,  the 
minimum-distance  estimate  is  almost  always  difficult  to  find 
(29:287-288) . 

The  Gram-Char lier  type  A  series  is  sometimes  used  to 
approximate  the  p.d.f.  of  a  distribution  whose  range  is 
doubly  infinite  (i.e.  f(x)>0  for  -oo<x<»).  It  is  based  on 
the  normal  distribution  and  its  derivatives  and  uses  a 
series  expansion  involving  Hermite  polynomials  to  approxi¬ 
mate  the  unknown  p.d.f.  (37:257-262;  9:222-227;  17:46-60). 
The  Gram-Char lier  type  B  series  is  based  on  the  Poisson 
distribution  and  involves  Poisson-Charlier  polynomials 
(17:72-81)  . 

In  a  recent  dissertation.  Hill  (17)  suggested  several 
ways  to  estimate  a  p.d.f.  if  a  finite  number  of  moments  or 


the  moment  generating  function  is  known.  If  the  moment 
generating  function  is  known  and  the  function  is  continuous 
over  the  positive  real  line,  then  the  p.d.f.  can  be  found  by 
finding  the  inverse  Laplace  transform  of  the  moment  gen¬ 
erating  function  (17:92-106).  Alternatively,  the  moment 
generating  function  could  be  used  to  obtain  moments  of  the 
distribution  and  one  of  the  series  expansions,  using  the 
Gram-Charlier  type  A  or  B  series  or  Laguerre  series,  could 
be  used  to  approximate  the  p.d.f.  (17:46-81).  Hill  also 
suggested  using  the  moments  to  fit  a  curve  of  the  Pearson 
family  of  frequency  distributions  (17:82-92). 

Technique  Selection 

The  H-function,  defined  in  Chapter  2  as 


m  n 

H  [x: { (ai,Ai)} ,  i-1, 

p  q 


,p;  ( (bj ,B j ) } ,  j«l, 


has  2(p+q)  parameters  to  be  estimated.  For  statistical 
distributions,  two  additional  parameters  are  included  to 
allow  for  scaling  and  to  ensure  that  the  H-function  distrib¬ 
ution  integrates  to  one  over  the  appropriate  range.  Thus, 
for  H-function  distributions,  there  are  2(p+q+l)  parameters 
to  be  estimated.  We  therefore  needed  a  method  of  estimating 
parameters  that  could  produce  estimates  for  a  vector  of 
parameters. 


Method  of  Maximum  Likelihood. 


As  noted  in  a  previous  section,  the  method  of  maximum 
likelihood  can  be  used  with  a  vector  of  parameters  and  is 


widely  used  for  this  purpose.  The  maximum  likelihood  esti¬ 
mates  also  possess  many  desirable  properties  of  estimators. 
Therefore,  we  attempted  to  obtain  maximum  likelihood 
estimates  for  the  parameters  of  the  H-function. 

Let  x1,X2/.../Xr  be  a  random  sample  from  the  H-function 
distribution 


m  n 

kH  [ex: { (ai,Ai) } ,  i=l,...,p;  {(bj,Bj)},  j=l,...,q], 
P  q 


m  n 

hereafter  abbreviated  as  kH  (cx).  Our  objective  is  to 

p  q 

obtain  point  estimates  for  the  2(p+q+l)  parameters  k,c,a^ 
and  A±  ( i=l,  ...,p) ,  and  bj  and  Bj  (j  =  l,...,q)  using  the 
method  of  maximum  likelihood. 

The  likelihood  function  is  simply  the  product  of  the 
individual  densities.  For  a  random  sample  of  size  r, 


r  m  n 

L(k,c,a^,  A^  f  u.jiB^|x^,...i  Xj.)  =  [  i  kH  ( cx^ ) 

h=l  p  q 


This  function  must  be  differentiated  with  respect  to  each  of 
the  parameters  k,c,a^,A^fbj,Bj.  Differentiation  with 
respect  to  k  or  c,  while  not  trivial,  is  easy  when  compared 
to  differentiation  with  respect  to  the  other  parameters. 


Mathai  and  Saxena  (26:19)  give  results,  due  to 
Buschman,  for  the  Mellin  transform  of  the  partial  deriva¬ 
tives  of  an  H-function  with  respect  to  its  parameters. 
These  results  imply  that  the  partial  derivatives  may  be 
brought  through  the  contour  integral  and  evaluated  using  the 
chain  rule.  For  example, 


1  0 

dH0  i(x) 
§b 
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P(b+Bs)  ^(b+BsJx-3 


ds 


and 


1  0 

dH  (x) 
0  1 

§B 


1 

2^i 


s  P(b+Bs)  'Hb+Bs ) x“s  ds 


Consider  the  special  case  where  k=c=m=q=l  and  n=p=0. 
Then  the  likelihood  function  is: 


r  10 

L(b , B I x^, . • • ,  xr )  =  [  |  H  ( x^) 

h=l  0  1 


r  I 
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f  P (b+Bs )  (xh)"s  ds 

-  TT  | 
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The  maximum  likelihood  estimates  for  b  and  B  are  the 
solutions  to  the  two  equations: 
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2TTi  j  P(b+Bs)^(b+Bs)  (xh)”sds 

[  TT  2'n‘i  f  P(b+Bs  )  (  x  ■  )_s  ds]  j  = 
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(3.1) 


Y  2m  s  ["(b+Bs)  ^(b+Bs)  (xO“sds 
h=l  JC 


[  JT  2"-i  /  f(b+Bs)  (xi)“S  ds3  =  0 


(3.2) 


Contour  integrals  are  usually  evaluated  by  summing  the 
residues  at  the  poles  of  the  integrand.  For  the  first 
contour  integral  in  Eq  (3.1),  zero  and  all  the  negative 

integers  are  poles  of  order  two.  Therefore,  it  could  be 

oo 

replaced  by  an  infinite  sum  of  residues,  say  X!  gtx^bjB, J) . 

J=0 

The  other  contour  integrals  in  Eq  (3.1)  are  simply  H- 
functions  which  equal 


iu.)B 


For  the  first  contour  integral  in  Eq  (3.2),  zero  and  all  the 
negative  integers  are  poles  of  order  two.  Therefore,  it 
could  also  be  replaced  as  an  infinite  sum  of  residues,  say 


y  u(xh,b,B,J).  Thus,  Eq  s  (3.1)  and  (3.2)  could  b< 
J«0 

rewritten  as: 
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The  maximum  likelihood  estimates  are  the  values  of  b  and  B 
which  satisfy  the  above  equations. 

It  should  be  noted  that  the  convenient  assumptions 
k=sc=m=q=l  and  n=p=0  caused  several  simplifications.  First, 
there  are  only  two  equations  to  be  solved.  Second,  the 
integrand  in  each  contour  integral  was  relatively  simple. 
With  more  complicated  integrands,  the  evaluation  of  contour 
integrals  by  the  sum  of  residues  becomes  more  difficult. 
Finally,  in  the  more  general  case,  the  H-functions  in  Eqs 
(3.1)  and  (3.2)  could  not  be  expressed  in  the  closed  form: 


1 

b  B 

i(x  )B  e"(  51 
B  j 


Instead,  each  would  be  an  infinite  sum  of  residues,  say 

oo 

Y  l(x^,b,B,J).  This  would  make  Eq  (3.1)  of  the  form: 

T _ J 


CO 

/  00  \1 

Y  g  ( xh ,  b ,  B ,  j ) 

TT  (  Y  i  ( x-i  ,b,B,  J)  ) 

J=0 

.  j^h  \  J=0  J  / J 

Since  no  general  results  are  known  for  the  product  of 
contour  integrals  or  the  product  of  infinite  series,  more 
research  is  required  before  maximum  likelihood  estimates  for 
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the  parameters  of  the  H-function  can  be  developed.  Thus, 
the  method  of  maximum  likel ihood  cannot  yet  be  used  to 
obtain  estimates  for  the  H-function  parameters. 


Method  of  Moments. 

The  method  of  moments  seems  to  be  the  next  most  widely 
used  technique  for  density  estimation,  in  spite  of  the 
concerns  mentioned  earlier.  As  seen  in  Chapter  2,  the 
analytic  moments  of  the  H-function  are  simply  products  and 
quotients  of  gamma  functions,  where  the  argument  of  each 
gamma  function  is  a  linear  combination  of  a  pair  of  the 
parameters.  These  moments  could  be  set  equal  to  the  moments 
of  the  data,  creating  as  many  equations  as  there  are  param¬ 
eters  to  be  estimated.  These  equations  can,  in  theory,  be 
solved  to  obtain  estimates  for  the  parameters.  The  method 
of  moments,  as  applied  to  H-functions,  will  be  fully  derived 
in  Chapter  5. 

The  other  parametric  density  estimation  techniques 
generally  do  not  produce  estimates  as  good  as  those  from  the 
method  of  moments.  Although  the  method  of  linear  least 
squares  can  also  produce  good  estimates,  the  H-function 
cannot  be  transformed  into  an  equation  that  is  linear  in  its 
parameters . 


IV  Nonlinear  Solution  Technique 


The  method  of  moments  produces  a  system  of  nonlinear 
equations  that  each  involve  products  and  quotients  of  gamma 
functions.  We  need  to  find  a  solution  of  these  simultaneous 
equations.  Numerical  analysis  techniques  already  exist  that 
can  accomplish  this  task.  A  survey  of  the  most  common 
techniques  was  conducted  to  determine  which  might  best  be 
applied  to  our  specific  problem.  Three  of  the  main  consid¬ 
erations  were  convergence  conditions,  amount  of  calculation, 
and  rate  of  convergence. 

Fixed-Point  Iteration 

In  fixed-point  iteration,  we  consider  two  nonlinear 
equations  in  the  form 

f(x,y)  =  o 

(4.1) 

g(x,y)  =  0 

We  rewrite  these  equations  by  taking  an  x  out  of  the  first 
equation  and  a  y  out  of  the  second.  This  gives  the 
equivalent  form 

x  =  F( x, y ) 

(4.2) 

y  =  G(  x,  y) 

so  that  any  root  (x,y)  will  solve  both  sets  of  equations. 
It  should  be  noted  that  there  are  many  ways  to  rewrite 
equations  from  the  form  of  Eq  (4.1)  to  the  equivalent  form 
of  Eq  (4.2). 
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The  fixed-point  iteration  begins  with  an  initial  guess 
(x0#y0)  and  generates  successive  approximations  from  the 
recursive  relationship  [6:84] 


xi+l  =  F(xi'Yi) 

Yi+1  =  G<xi'Yi> 


(4.3) 


It  is  possible  to  accelerate  the  iteration  process  by  using 
the  most  current  information  on  in  the  second  equation 

of  Eq  (4.3).  This  produces 

xi+l  =  F(xi^i) 

Yi+1  =  G^xi+l'Yi) 


which  will  converge  provided  the  original  iterative  process 
in  Eq  (4.3)  converges  [5:446]. 

For  the  original  iterative  process,  convergence  occurs 
under  the  following  sufficient  (but  not  necessary) 
conditions  : 

1.  F  and  G  along  with  their  first  partial  derivatives 
are  continuous  in  a  neighborhood  about  the  root 

(x.y) , 


2.  For  all  points  in  the  neighborhood. 


IdF 
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aci 

lax 
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ay 
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<L  M  <  1 

£  M  <  1 


for  some  M,  and 

3.  The  initial  approximation  (Xp,yQ)  is  taken  from 
the  same  neighborhood  [6:84;  39:130;  35:223]. 
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Moreover,  if  M  is  very  small  for  all  steps  in  the  iteration 
(Xi/Yi),  then  the  iteration  converges  quickly  relative  to 
the  case  when  the  magnitude  of  M  is  near  one  [39:130]. 

We  define  the  order  of  convergence  as  a  measure  of  the 
speed  or  rate  of  convergence.  Order  of  convergence  is  the 
lowest  value  of  n  such  that  the  n^-derivative  of  g(x) 
evaluated  at  the  solution  x  does  not  equal  zero.  For  this 
reason,  fixed-point  iteration  of  the  type  x  =  g(x) 
generally  has  first-order  or  linear  convergence  [39:84-85]. 
Newton's  Method 

Next,  we  discuss  another  fixed-point  iteration  of  the 
type  x  =  q( x )  which  has  second-order  or  quadratic  conver¬ 
gence  [44:145-146].  In  this  case,  g(x)  is  chosen  so  that 
its  first  derivative  vanishes  at  the  solution  x. 

The  formula  for  a  Newton  iteration  can  be  derived 
through  a  Taylor  series  expansion  of  f(x).  When  f(x)  is 
twice  continuously  differentiable,  then 


\  et  \  f'U.)U.  -x.)  ^  f '  '  (  £  )  (x  -x .  ) 2 

f(xi+i)  =  f  (x±)  +  _ 1  1+1  1  +  1+1  1 


11 
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where  £  lies  between  x^  and  x^  +  1.  Suppose  that  x^+1  is 
chosen  so  that  f(x^+^)  is  nearly  equal  to  zero  and  also  that 
(x^  +  1-Xj_)2  is  sufficiently  small  so  that  the  last  term  can 
be  neglected.  Then  the  above  Taylor  series  simplifies  to 


0  =  f  (X  •  )  +  f  '  (X-  )  (X  •  —X;  ) 


relationship  for  Newton's  Method  as  [44:132] 


x 


i+1 


f(x±) 


(4.4) 


For  the  general  form  x  =  g(x),  the  iterative  process 
in  Eq  (4.4)  will  have  second-order  convergence  if  the 
derivative  of  g(x)  is  zero  at  the  solution  x.  By- 
differentiating, 

g'(x)  =  f(x)  f"(x) 

[  f 1  ( x )  ]  2 

and  at  x  =  x, 

g'(x)  =  f(x)  f"(x) 

[ f ' ( x) ]2 

Since  f(x)=0,  the  numerator  is  equivalent  to  zero.  If  the 
root  x  is  simple,  then  f'(x)  is  nonzero,  g'(x)  is  zero,  and 
the  iteration  converges  with  order  two.  However,  if  multi¬ 
ple  roots  exist  at  x,  then  L'Hospital's  Rule  must  be  used  to 
show  the  convergence  degenerates  to  first-order  [39:85-89; 
44:132-133]. 

If  the  derivative  of  f(x)  cannot  be  solved  explicitly, 
then  it  can  be  solved  using  the  approximation  formula 

f'(x)  ~  £  L  2£  .+.  A .  1 — ~  ,^(  x)  (4.5) 

€ 

for  very  small  5  [6:278]. 


To  find  the  solution  to  the  system  of  n  nonlinear 


equations 

fx  (St)  =  0 
f2(St)  =  O 

(4.6) 

fn(X)  -  o 

the  recursive  relationship  in  Eq  (4.4)  can  be  expanded  to 


\ 


*i+1  -  -  J(2i)"1  F(X±) 


(4.7) 


where  J(X^)  is  called  the  Jacobian  matrix  defined  by 


J(S)  = 


3fl<*> 

3*1 

3£2(X) 

a*! 


df  (X) 
c*x2 

af2(x) 

dx9 


df X(x) 


n 


df2(f> 

dxn 


df  (X)  df  (X) 
n  n 

dxi  dx2 

evaluated  at  the  ith  iteration 
the  equations  in  Eq  (4.6)  at  X^  is  F(X^)  [5:450]. 

The  approximation  formula  for  derivatives  in  Eq  (4.5) 
can  also  be  expanded  to  become 


df  (X) 
n 


of  X  and  the  evaluation  of 


df,(x4)  „  f,<x,  +  -  f,(x.) 


(4.8) 


where  £  is  very  small  and  e^  is  a  n-vector  whose  only  non¬ 
zero  entry  is  a  one  in  the  ktn  row  [5:456;  30:26]. 

For  the  recursive  relationship  in  Eq  (4.7),  convergence 
occurs  under  the  following  sufficient  (but  not  necessary) 
conditions  : 

1.  f l » f 2'  *  ^n  a^on9  with  all  derivatives  through 
second  order  are  continuous  in  a  neighborhood 
about  the  root  vector  x, 

2.  The  Jacobian  J(2^)  does  not  vanish  in  the  same 
neighborhood,  and 

3.  The  initial  approximation  StQ  is  chosen  suffi¬ 
ciently  close  to  the  root  vector  x  [6:86;  5:449]. 

When  the  generalized  Newton's  method  converges,  it  has 
second-order  convergence  if  the  roots  are  simple.  But  it  is 
difficult  to  insure  that  the  determinant  of  the  Jacobian  is 
not  zero.  Furthermore,  the  necessity  to  invert  the  Jacobian 
at  each  iteration  requires  many  computations  and  thus 
simpler  methods  are  needed  in  most  cases  [39:133]. 

Modified  Newton's  Method 

For  two  nonlinear  equations  in  two  unknowns,  the  modi¬ 
fication  of  Newton's  method  consists  of  applying  the  single¬ 
variable  Newton  method  two  times,  once  for  each  variable. 
Each  time  this  is  done,  the  other  variable  is  assumed  to  be 
fixed.  Succeeding  approximations  are  then  generated  from 
the  recursive  relationship  [39:136] 


‘i+1 


=  x. 


-  f(xi.,yi) 


V: 


fxU“Tyi) 

g(xi,yi) 


(4.9) 


which  can  be  accelerated  by  using  the  most  current 
information  on  in  the  second  equation  of  Eq  (4.9) 

[6:88] . 

Note  that  generally  we  could  use  either  f  or  g  to 
calculate  the  new  x  and  use  the  other  function  to  calculate 
y.  One  of  these  choices  will  usually  converge  while  the 
other  diverges,  depending  on  the  exact  problem.  For  n 
nonlinear  functions  in  n  unknowns,  there  are  nl  ways  of 
choosing  the  n  functions  to  find  n  unknowns.  Often  only  one 
of  these  choices  will  converge  and  this  is  the  main  disad¬ 
vantage  of  this  method  [39:136-143].  But  if  the  correct 
combination  is  found,  the  convergence  rate  will  be  remark¬ 
ably  rapid  and  faster  than  linear  convergence  [6:87-89]. 
However,  we  need  to  find  another  simple  method  with  a  better 
chance  for  convergence  to  the  root. 

Steepest  Descent  Method 

Next,  we  discuss  a  gradient  search  technique.  In 
steepest  descent,  the  recursive  relationship  for  a  general 
system  of  nonlinear  equations  is 

*i+1  -  -  J F (X±) 

which  is  similar  to  Eq  (4.7)  except  that  the  transpose  of 
J( St )  is  used  in  place  of  the  inverse  [32:62-63].  Another 
way  to  look  at  steepest  descent  is  to  consider  the  problem 
of  minimizing 
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M(*jc>  =  I  Cfi(Sk)]2 

i=l 


(4.10) 


The  value  of  2k  which  causes  M(Xk)  in  Eq  (4.10)  to  equal 
zero  will  also  solve  the  original  set  of  n  nonlinear  equa¬ 
tions  [24:244-245].  Although  it  only  promises  linear  con¬ 
vergence,  "the  method  of  Steepest  Descent  has  been  found  to 
be  an  effective  way  of  getting  reasonably  close  to  the  solu¬ 
tion"  [32:63].  For  this  reason,  hybrid  algorithms  based 
initially  on  the  method  of  Steepest  Descent  followed  by 
Newton's  method  can  be  very  reliable  for  systems  of  non¬ 
linear  equations  [14:36-37]. 

Levenberg-Marguardt  Method 

Such  hybrid  algorithms  are  quite  commonly  used  today. 
When  dealing  with  a  system  of  n  nonlinear  equations,  the 
recursive  relationship  for  the  Levenberg-Marquardt  method  is 
given  by 

*i+1  “  Xi  "  H<*i>  F<*i> 

where  H(2i)  -  [ J(2i)TJ(5?i)  +  AiI]~1J(Xi)T  and  F(2i)  is  the 
evaluation  of  the  equations  in  Eq  (4.6)  at  5?^ 
[32:63].  Note  that  as  A^  increases,  the  step  vector  H(X^) 
tends  toward  the  pure  Steepest  Descent  vector.  On  the 
other  hand,  if  A^=0,  then  the  method  reduces  to  Newton's 
method.  By  reducing  A^  systematically,  the  hybrid  itera¬ 
tion  combines  the  better  features  of  both  methods  [32:63-64] 
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However,  the  problem  of  inverting  and  solving  the 
Jacobian  matrix  at  each  iteration  is  still  a  major  weakness 
because  the  Levenberg-Marquardt  method  still  involves  an 
inversion.  The  number  of  computations  is  excessive  even  for 
a  computer.  The  uolution  to  these  problems  involves  algo¬ 
rithms  that  are  known  as  quasi-Newton  [4:577-578].  The  idea 
of  a  quasi-Newton  algorithm  is  to  eliminate  the  calculations 
involved  with  the  inversion  of  the  Jacobian  matrix  [14: 
38-45;  32:577-578]. 

Broyden's  Procedure 

Broyden's  procedure  can  be  used  in  an  iteration  method 
to  avoid  the  inversion  of  the  Jacobian  matrix  at  every 
single  iteration.  The  approximation  to  the  inverse  of  the 
Jacobian  matrix  A(£)  is  updated  at  each  iteration  using  the 
formula 


A(Xi+1)  =  A(Sti)  + 


C(Si+i  -  A(Xi)Yi4.1)Si+1  A(3E±)  ] 
Sill  M*i>  Yi+1 


where  Yi+1  =  F<Xi+i)  ”  F(*i>  and 

-  2^  [5:455-460;  4:581]. 

Of  course,  for  the  first  iteration  the  actual  Jacobian 
matrix  J(X0)  must  be  found  explicitly  or  approximated  using 
Eq  (4.8).  Tnen  J(ZQ)  must  be  inverted  once  before 


Broyden's  procedure  can  begin.  After  that,  the  procedure 


always  produces  an  approximation  to  the  inverse  of  the 


Jacobian  matrix  A(X^)  which  can  be  used  to  replace 
J(2^)-^,  the  actual  inverse  of  the  Jacobian  matrix.  The 
recursive  relationship  is 

*i+l  "  *i  "  A(X±)  F<Xi> 

which  is  the  same  as  Eq  (4.7). 

This  procedure  significantly  reduces  the  number  of 
arithmetic  calculations  and  still  provides  super  linear 
convergence  [16:5-6].  Therefore,  iterative  methods  with 
second-order  convergence  will  approach  second-order  conver¬ 
gence  when  Broyden's  procedure  is  used.  Quadratic  conver¬ 
gence  will  be  obtained  as  the  approximation  to  the  inverse 
of  the  Jacobian  matrix  A(X^)  becomes  better  [5:456], 

Powell 1 s  Method 

The  best  combination  of  methods  and  procedures  studied 
so  far  would  be  a  quasi-Newton  hybrid  algorithm.  This  was 
found  in  an  IMSL  routine  named  ZSPOW  which  contains  M.J.D. 
Powell's  hybrid  method  (HYBRD1 )  for  nonlinear  equations 
[32:87-114].  ZSPOW  not  only  includes  the  beneficial  fea¬ 
tures  of  the  Levenberg-Marquardt  method,  but  also  implements 
the  calculation-saving  strategy  of  Broyden's  procedure 
[16:6-7;  30:45]. 

In  a  comparison  of  available  software  which  solves 
systems  of  nonlinear  equations,  HYBRD1  had  outstanding  per¬ 
formance.  Also,  initial  estimates  of  the  parameters  had 


little  effect  on  the  convergence  [16:24,  41-44].  All  of 
this  comes  while  still  providing  super  linear  convergence. 
Therefore,  the  method  used  in  ZSPOW  will  nearly  obtain 
second-order  convergence. 

In  summary,  the  IMSL  routine  named  ZSPOW  has  several 
advantages.  First,  ZSPOW  is  a  hybrid  routine  which  permits 
a  bad  initial  guess  of  the  root.  Second,  ZSPOW  is  quasi- 
Newton  with  a  convergence  rate  that  is  nearly  second-order. 
Next,  the  user  needs  to  supply  only  a  subroutine  that  con¬ 
tains  the  system  of  nonlinear  equations  that  has  to  be 
solved.  Finally,  ZSPOW  outputs  error  messages  if  the 
iteration  does  not  make  good  progress. 

We  applied  this  numerical  analysis  technique  to  our 
system  of  nonlinear  equations.  The  program  described  in  the 
next  chapter  uses  ZSPOW  to  produce  accurate  estimates  for  up 
to  ten  unknown  variables  from  the  same  number  of  equations. 


V  Methodology 

The  method  of  moments  has  been  selected  as  the  most 
appropriate  curve- fitting  technique  to  estimate  the  param¬ 
eters  of  the  H-function.  The  method  involves  equating  the 
appropriate  number  of  analytic  H-function  moments  with  the 
same  number  of  data  moments.  This  will  create  as  many 
equations  as  there  are  parameters  to  be  estimated.  As  seen 
in  the  following  derivation,  the  number  of  equations  can  be 
reduced  by  two  through  algebraic  manipulation. 

Generation  of  Equations 

As  discussed  earlier  in  Chapter  2,  the  rtlri  moment,  /jr, 

m  n 

of  the  H-function,  kH  (cx),  is  defined  by  the  following 

P  q 

equations 

^r+T  Kr+l) 

where 


I(r+1) 


m 

TT 

3=1 

P 

TT 

i=n+l 


r ( bj+B j+B jr )  JT  TU-ai-Ai-Air) 

_  q  “ 

['{*i+&i+Aix)  T T  rd-b1-Bj-B1r) 

j=m+l  J 


Since  I(r+1)  has  two  parameters  in  the  argument  of  every 
gamma  function,  each  H-function  moment  equation  will  involve 
2(p+q)+2  unknown  parameters.  The  same  number  of  equations 


needs  to  be  generated  by  setting  |ir  equal  to  the  data's  rth 
moment . 

Since  consecutive  data  moments  should  be  used,  the 
variable  r  could  take  on  the  values  from  1  to  2(p+q)  +  2.  Let 
the  r*’*1  moment  of  the  data  be  represented  by  Mr.  By  solving 
each  equation  for  k, 


k  =  Mr  c 

iTr+rr 

for  r=l,  2, 2(p+q)+2.  Since  all  the  equations  equal  k. 


k  =  MX  c2  .  M2  =3  .  _  .  «2(_p+q)-f2  c2<P+q)«  (5-1) 


iXTi~  JXTT 


I ( 2 (p+q)+3 , 


The  adjacent  equations  in  Eq  (5.1)  can  be  solved  for  c  to 
give 


c  „  Mr  I(r+2) 

M  ,  I ( r+1 ) 
r+1 

for  r=l,  2, ...,  2(p+q)  +  l.  Since  all  the  equations  equal  c. 


c  =  M1  =  M2  i(4)  = 

M  I  ( 2 )  M""lf3) 

2  3 


...  =  M2(p+q)+l  I(2(ptq)-*-3)  (5.2) 


M  ,  „  I (2 ( p+q) +2 , 
2(p+q)+2 


The  adjacent  equations  in  Eq  (5.2)  can  be  solved  to  give  the 
following  homogeneous  equations: 


MiMi+2  Cl(i+2)]2 
(Mi+1)2  I ( i+1 )I ( i+3 ) 


-1  =  0 


(5.3) 


for  i=l , 2, . . . , 2 (p+q) . 
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Note  that  although  this  algebraic  manipulation  has 
reduced  the  number  of  equations  from  2(p+q)+2  to  2(p+q),  the 
analyst  will  still  be  required  to  calculate  all  2(p+q)+2 
data  moments.  Note  also  that  if  the  zeroth  moment  is  used, 
then  Eq  (5.3)  will  be  evaluated  for  i=0, 1, ...,  2(p+q)-l. 

An  example  of  this  equation  generation  technique  may  be 
helpful  at  this  time.  For  the  generalized  gamma 
distribution,  the  H-function  is  given  as 

1  O 

kh  [cx : (b,B) ] 

0  1 

Since  p+q=l,  the  method  should  generate  four  equations  and 
then  reduce  the  number  of  equations  by  two.  Since  l(r+l)= 
P(b+B(r+1)),  the  four  equations  produced  by  the  method  of 
moments  will  be 

M1  =  Hb+ZB) 

M2  =  T (b+3B) 

M3  =  r(to+4B) 

c 

M4  =  ^P(b+5B) 


By  solving  for  k  and  setting  the  equations  equal  to  each 
other. 


p( b+2B  J  P(b+3B)  r( b+4B)  r( b+5B) 


The  adjacent  equations  are  then  solved  for  c  and  set  equal 
to  each  other  to  produce 

c  =  M1  Rb+3B)  _  M2  f(b+4B)  =  M3  P ( b+  5 B ) 


The  adjacent  equations  are  again  solved  and  the  result  is 
the  two  following  homogeneous  equations: 


M1M3  [ p( b+3B ) ] 2 


(M2)2  p(b+2B)P(b+4B) 
M2M4  [P(b+4B)]2 


(M3)  P(b+3B)P(b+5B) 


-1  =  0 


Once  the  two  equations  have  been  solved  for  the  parameters  b 
and  B,  then  the  estimates  may  be  used  in  Eq  (5.5)  to  solve 
for  the  value  of  c.  After  this  has  been  done,  k  can  be 
found  using  Eq  (5.4). 

As  expected  in  the  example,  the  method  generated  two 
equations  in  two  variables.  The  two  equations  still  require 
four  data  moments  to  be  evaluated.  To  emphasize  a  point 
previously  made,  if  MQ  through  M3  had  been  used  instead, 
then  the  method  would  have  resulted  in  the  different 
homogeneous  equations: 


MqM2  CP(  b+2B ) ] 2  _  x 

(Mx)2  p(b+B )P( b+3B) 

MxM3  CP(b+3B)]2  _  1 

(M2)2  P*(b+2B)P(b+4B) 


(5.6) 


These  two  equations  are  in  the  form  of  Eq  (5.3)  when  i=0  and 
when  i=l  respectively. 

Program  Development 

The  general  form  of  the  equations  in  Eq  (5.3)  verifies 
the  fact  that  a  system  of  nonlinear  equations  needs  to  be 
solved.  Even  the  simple  case  demonstrated  in  the  example 
produced  two  simultaneous  nonlinear  equations.  Chapter  4 
concluded  that  the  IMSL  routine  named  ZSPOW  would  be  the 
best  numerical  analysis  technique  to  use  for  nonlinear  sys¬ 
tems  of  equations.  Therefore,  our  program  implements  ZSPOW 
in  order  to  find  estimates  for  the  desired  parameters. 

Phase  One . 

The  actual  program  development  was  modularized.  That 
is,  the  program  was  divided  into  four  successive  phases  so 
that  each  was  easier  to  solve.  The  first  phase  involved 
programs  that  could  take  perfect  moments  from  known  distrib¬ 
utions  and  use  ZSPOW  to  obtain  estimates  for  the  parameters. 
Only  the  aj_,  A^,  bj,  and  Bj  need  to  be  guessed  initially. 
The  values  of  c  and  k  can  be  found  using  the  formulas  in  Eq 
(5.2)  and  Eq  (5.1)  respectively. 

The  early  programs  in  this  phase  could  only  handle  one 
special  case  of  the  H-function  at  a  time.  This  was  because 
ZSPOW  required  a  subroutine  to  contain  the  equations  which 
needed  to  be  solved  and  the  early  programs  only  used  the 
actual  equations  like  those  in  Eq  (5.6).  In  addition,  this 
phase's  early  programs  could  only  output  the  results  for  the 


one  special  H-function.  However,  these  early  runs  did  dem¬ 
onstrate  the  merit  of  the  ZSPOW  numerical  analysis  technique 
for  finding  accurate  estimates  of  the  parameters. 

The  later  runs  in  phase  one  became  more  general.  These 
programs  could  generate  equations  in  the  form  of  Eq  (5.3) 
and  output  the  results  for  the  input  values  of  m,  n,  p,  and 
q.  The  only  restriction  was  that  p+q  had  to  be  smaller  than 
six  due  to  matrix  dimension  limitations.  Our  research  was 
not  expected  to  go  past  a  2nd  order  H-function.  Also,  the 
input  variable  FLAG  controlled  whether  the  zeroth  moment  was 
used.  Therefore,  the  user  could  generate  the  correct 
equations  by  using  the  proper  input  of  FLAG. 

This  concluded  the  first  phase  programs.  Accurate 
results  were  achieved  with  perfect  moments  for  all  1st  and 
2nd  order  H-functions.  At  times  when  the  initial  guess  of 
the  parameters  was  far  away  from  the  proper  number,  ZSPOW 
would  not  converge  on  the  expected  root.  This  problem  was 
corrected  during  the  second  phase. 

Phase  Two. 

The  second  phase  involved  attempts  to  control  the  ini¬ 
tial  guess  of  the  parameters.  The  first  control  consisted 
of  checking  the  initial  guess  against  a  set  of  H-function 
convergence  conditions  (26:3;  7:72).  Some  of  these  condi¬ 
tions  were  mentioned  previously  in  Chapter  2,  but  now  we 
fully  discuss  the  subject. 
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If  a^,  bj,  ,  Bj,  m,  n,  p,  and  q  hold  their  usual 
meaning  in  the  definition  of  the  H-function,  then  the  new 
terms  D,  E,  and  L  can  be  defined  by  the  following  equations: 


Chapter  2  also  pointed  out  that  all  Left  Half-Plane  (LHP) 
m 

poles  of  TT  n^+B-jS)  must  lie  to  the  left  of  C,  and  all 
j=l  J  J 

n 

Right  Half-Plane  (RHP)  poles  of  XT  Pd-a-A.^)  must  lie  to 

i=l 

the  right  of  C^.  From  this  point  on,  we  refer  to  the  path 
of  integration  ( C )  as  the  w  line.  Since  there  may  be  a 
significant  distance  between  the  right-most  LHP  pole  and  the 
left-most  RHP  pole,  the  u>  line  may  be  placed  anywhere 
within  the  range  defined  by  those  two  values.  This  distance 
will  be  referred  to  as  the  w  range  where  Wi0w  '*'s  the 
right-most  LHP  pole  and  whigh  is  the  left-most  RHP  pole. 

With  these  definitions,  a  completely  specified  H- 
function  represented  by  an  infinite  sum  of  residues  does  not 
converge  under  any  of  the  following  conditions: 


Case  1.  D  =  0  E<0  L>(E)  a>low 

Case  2.  D  =  0  E>0  L>(E)  w^g^ 

Case  3.  D  =  O  E=0  L>0 

Case  4.  D  =  0  E  =  0  L=0 

Case  5.  D  <  O 

Case  6.  «low  >  «high 

For  all  other  cases,  the  H-function  will  converge. 

This  convergence  check  was  placed  in  a  program  sub¬ 
routine.  After  the  program  user  makes  an  initial  guess  of 
the  parameters,  the  completely  specified  H-function  is 
checked  against  the  convergence  conditions.  This  forces  the 
initial  guess  to  meet  convergence  conditions  and  helps  ZSPOW 
find  the  correct  root.  As  a  last  check,  the  final  estimate 
of  the  parameters  determined  by  ZSPOW  is  also  checked  for 
convergence . 

Up  to  this  point,  the  user  had  been  given  the  responsi¬ 
bility  for  the  initial  guess  of  the  parameters.  The  second 
control  on  the  initial  guess  involved  supplying  initial 
guess  default  values.  Two  requirements  were  built  in. 
First,  the  default  values  were  designed  to  always  meet  the 
initial  convergence  condition  check.  Second,  the  default 
guesses  only  take  on  values  between  zero  and  one.  The 
reason  for  this  may  not  be  clear.  For  common  distributions, 
when  an  H-function  parameter  is  specified  by  an  exact  number 
that  number  will  take  on  the  value  0,  1/2,  or  1  (Table  I). 
When  a  H-function  parameter  is  represented  by  a  variable,  no 


■-  ^  -  %  % 


default  guess  will  be  appropriate  for  all  possible  values  of 
the  variable.  Therefore,  the  default  guess  of  the 
parameters  is  as  good  as  possible. 

Many  runs  on  1st  and  2n<^  order  H-functions  were  per¬ 
formed  at  this  time  to  validate  the  use  of  the  variables 
FLAG  and  GUESS.  When  FLAG  =  0,  the  zeroth  moment  was  used. 
When  FLAG  =  1,  the  moments  began  with  the  first  moment.  The 
same  results  were  achieved  on  each  run  no  matter  how  FLAG 
was  set.  Since  we  were  concerned  with  the  inaccuracy  of 
higher  degree  moments,  most  future  runs  were  performed  with 
FLAG  =  0.  This  is  the  suggested  configuration  for  running 
the  program.  However,  when  the  zeroth  moment  is  inaccurate, 
the  program  user  may  wish  to  set  FLAG  =  1. 

If  GUESS  =  0,  the  user  supplied  the  initial  guess  of 
the  H-function  parameters.  When  GUESS  =  1,  the  default 
initial  guess  was  used.  The  default  guess  also  performed  as 
expected.  If  distributions  had  actual  parameter  values 
around  the  range  (0,1),  then  GUESS  =  1  converged  to  the 
correct  root.  If  distributions  were  run  that  had  actual 
parameter  values  too  far  from  the  (0,1)  range,  GUESS  =  1 
produced  an  error  message.  For  example.  Beta  (9=2,0=10) 
will  not  run  using  the  default  guess  because  the  initial 
estimate  of  small  a=0.7  is  not  close  enough  to  the  actual 
value  of  a=11.0.  Therefore,  we  suggest  that  the  program 
user  set  GUESS  =  0  and  make  an  initial  guess  unless  he 


thinks  the  unknown  distribution  has  actual  values  for  the 


parameters  near  the  range  (0,1).  Table  II  lists  the  dis¬ 
tributions  which  may  always  be  run  with  the  variable  GUESS 
set  to  one. 


Table  II.  Distributions  With  Constant  H- function  Parameters 


(a,  A) 

(b,B) 

Exponential 

(0,1) 

Rayleigh 

_ 

(1/2, 1/2) 

Maxwell 

- 

(1,1/2) 

Half-Normal 

- 

(0,1/2) 

Uniform 

(1,1) 

(0,1) 

Bessel 

- 

(0,1/2)  (0,1/2) 

All  other  runs  with  GUESS  =  0  converged  when  a  reasonable 
guess  of  the  parameters  was  made. 

In  conclusion,  the  accurate  convergence  to  the  root  by 
ZSPOW  was  enhanced  by  the  convergence  checker  subroutine  and 
the  option  for  an  initial  default  guess  of  the  parameters. 
At  this  time,  we  had  produced  a  program  which  could  consis¬ 
tently  fit  more  special  functions  than  any  previous  general 
special  function  procedure. 


Phase  Three. 


In  all  previous  computer  runs,  exact  moments  had  been 
used.  In  this  respect,  the  computer  program  had  been  veri¬ 
fied.  By  verification,  we  mean  the  program  performed  as  it 
was  expected  to  perform.  Phase  4  would  continue  the  verifi¬ 
cation  process  by  insuring  that  error  messages  performed  as 
desired.  On  the  other  hand.  Phase  3  concentrated  on  valida¬ 
tion.  If  the  program  was  to  be  valid,  it  should  not  require 
perfect  moments. 

One  way  to  input  imperfect  moments  would  be  to  add  some 
error  to  each  moment.  A  more  reasonable  approach  might  be 
to  calculate  the  moments  from  raw  data.  But  an  analyst 
cannot  always  control  the  type  of  data  he  will  receive. 
Therefore,  Phase  3  added  the  capability  to  input  the 
following  four  types  of  data: 

Type  0.  Previously  calculated  moments 

Type  1.  Univariate  deviates 

Type  2.  Ordered  pairs  from  a  relative  frequency 

Type  3.  Ordered  pairs  from  a  continuous  function 

The  last  three  types  were  new  and  each  needed  a  subroutine 
that  could  calculate  moments  from  raw  data. 

For  type  one,  the  data's  r1*^1  moment  was  calculated  by 

n  (x.)r 


using 


where  there  are  n  univariate  data  points.  For  type  two,  the 


data's  r^  moment  was  calculated  by  using 


m 


M  =  Y  (X,)r  Pr(X=xi) 
r  j=l  11  3 


(5.8) 


m 


where  there  are  m  ordered  pairs  and  Y  Pr(X=x^)  =  1.  For 

3=1  J 


.th 


type  three,  the  data's  rttl  moment  is  exactly  represented  by 


M 


r  "  J  <x>r 
j-  00 


f(x)  dx 


This  formula  may  be  approximated  by  using 


m 


Mr  “  L  (xk^r  f(xk^  Ax 
k=l 


(5.9) 


m 

where  there  are  m  ordered  pairs,  Y  f(xv)  need  not  equal 

k=l 

one,  and  Ax  is  the  interval  between  the  xk  values.  Since 
Eqs  (5.8)  and  (5.9)  were  similar  and  both  involved  ordered 
pairs,  they  were  combined  into  a  single  subroutine. 

Recall  that  the  H-function  is  only  applicable  for  con¬ 
tinuous  functions.  This  means  that  even  though  discrete 
distributions  like  the  Poisson  can  generate  data  in  the  form 
of  type  one  or  two,  the  H-function  will  not  be  useful  in 
these  instances.  Only  univariate  data  from  continuous  dis¬ 
tributions  defined  over  positive  x  should  be  used  in  type 


one.  Similarly,  only  relative  frequency  data  from  contin¬ 
uous  distributions  defined  over  positive  x  should  be  used  in 
type  two. 

At  this  stage,  we  generated  data  to  validate  the  pro¬ 
gram.  Type  one  data  consisted  of  deviates  generated  from 
standard  distributions.  For  example,  the  exponential 
deviates  with  0=2  were  generated  by  IMSL  using  the  formula 

x  =  -[ In ( 1-z) 3/2 

where  z  is  a  random  number  from  a  uniform  (0,1) 
distribution . 

Type  two  data  was  created  by  classifying  the  type  one 
deviates  into  intervals.  This  produced  ordered  pairs 
(x,f(x))  where  x  is  the  midpoint  of  the  interval  and  f(x)  = 
Pr(x=x).  The  Pr(X=x)  is  the  proportion  of  deviates  within 
the  interval.  In  the  exponential  example,  the  ordered  pairs 
may  look  something  like 

( .5, .5) 

(1.5, .25) 

(2.5, .10) 

(3.5, .06) 

(4.5, .04) 

(5.5, .03) 

(6.5, .02) 

The  first  ordered  pair  represents  the  fact  that  50%  of 
the  exponential  deviates  were  located  within  the  interval 
(0,1).  In  truth,  an  analyst  would  never  combine  type  one 
data  into  intervals  to  create  type  two  data.  By  doing  this. 


he  would  only  lose  information.  We  go  from  type  one  to  type 
two  data  only  for  the  ease  of  generating  ordered  pairs  of  a 
relative  frequency  from  a  continuous  distribution. 

Finally,  category  three  data  was  created  by  assuming 
the  function  was  known.  This  produced  ordered  pairs 
(x,f(x))  which  basically  plot  the  function  at  a  limited 
number  of  points.  To  keep  the  mathematical  calculations 
simple,  the  function  is  evaluated  so  that  Ax  is  the  same 
between  all  x  values.  In  the  exponential  example, 

f ( x )  =  2e”2x 

and  the  ordered  pairs  would  be  created  by  observing  the 
value  of  f(x)  as  x  is  incremented  from  0.01  to  10  by  0.01. 

We  were  careful  that  each  type  of  data  was  generated 
from  the  same  distributions  and  functions  that  had  been  used 
for  verification.  The  moments  calculated  from  the  data  by  a 
moment  generating  program  were  compared  to  the  perfect  mo¬ 
ments.  The  data  moments  were  fairly  close  to  the  true 
moments  but  they  became  increasingly  inaccurate  for  higher 
degree  moments.  Next,  we  validated  the  convergence  of  the 
phase  three  program  with  the  inaccurate  moments.  The  pro¬ 
gram  converged  to  nearly  the  same  root  that  had  been  ob¬ 
served  during  verification.  Finally,  the  actual  raw  data 
was  input  into  the  program.  The  program  always  converged  to 
the  same  root  as  had  been  seen  with  the  inaccurate  moments. 
This  validated  our  program's  ability  to  handle  raw  data. 
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In  conclusion,  the  program  now  had  the  following  fea¬ 
tures.  The  program  could  input  either  previously  calculated 
moments  or  three  types  of  raw  data.  If  raw  data  was  input, 
the  program  would  generate  the  appropriate  number  of  data 
moments.  The  program  would  then  supply  a  default  initial 
guess  of  the  H-function  parameters  if  the  user  did  not 
supply  his  own.  The  completely  specified  H-function  would 
next  be  checked  against  the  convergence  conditions.  If  they 
were  not  satisfied,  the  program  would  be  terminated.  If  the 
convergence  conditions  were  met,  the  program  would  continue 
by  using  ZSPOW  to  get  final  estimates  of  the  H-function 
parameters.  Once  ZSPOW  had  finished,  the  final  estimates 
were  again  tested  against  the  convergence  conditions.  The 
final  results  were  then  output  along  with  an  error  message 
if  the  last  convergence  test  had  failed.  The  program's 
flowchart  can  be  seen  in  Figure  1. 

The  reason  for  the  output  even  if  the  final  estimate 
did  not  meet  convergence  conditions  was  to  check  if  the 
inaccuracy  of  the  moments  had  caused  the  problem.  An  exam¬ 
ple  of  this  might  be  an  estimated  B=1.0  in  the  numerator  and 
an  estimated  A=1.01  in  the  denominator.  This  would  produce 
a  D=*-.01  in  Eq  (5.7)  and  therefore  the  H-function  would  not 
meet  the  convergence  conditions.  However,  both  estimates 
should  probably  equal  one.  The  estimate  of  A  may  be  off  due 
to  imperfect  moments.  The  only  way  to  find  this  type  of 
error  is  to  always  output  any  final  results  of  ZSPOW. 
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'Moments 

or 

Raw  Data i 


Moments 

are 

Calculated 


/  User  Guess 
I  or 

Default  Guess 


However,  many  other  types  jf  err.jrs  ^ill  cause  the  program 

to  stop. 

Phase  Four. 

The  final  phase  added  error  checks  and  comments  to  the 
program.  The  comments  were  intended  to  help  the  user  under¬ 
stand  the  program.  The  error  messages  were  intended  to  keep 
the  user  within  the  bounds  of  the  program  and  to  give  the 
user  some  idea  of  why  the  program  was  stopped.  In  order  to 
get  the  most  information  from  the  error  messages,  it  is 
important  to  know  the  order  of  error  checks.  For  instance, 
if  three  conditions  are  checked  and  the  program  stops  be¬ 
cause  the  third  condition  is  not  satisfied,  then  the  user 
can  reason  that  the  other  two  conditions  were  satisfied. 
The  following  list  contains  the  successive  order  of 
necessary  conditions  for  the  program  to  continue  running: 

1.  FLAG  =  0  or  1 

2 .  0  <  m  <  q 

3 .  0  i  n  <  p 

4.  p+q  <  5 

5.  TYPE  =0,  1,  2  or  3 

6.  If  TYPE  =  0,  then  number  of  moments  =  2(p+q)+2 

7.  If  TYPE  =1,  2  or  3,  then  amount  of  data  >  20 

8.  If  TYPE  =  2  or  3  then  Ax  the  same  between  all  pairs 

9.  GUESS  =  0  or  1 

10.  Initial  guess  meets  convergence  conditions 


11. 


ZSPOW  runs  without  IMSL  errors 


In  addition  to  these  error  messages,  the  program  will  output 
an  error  message  if  the  final  estimate  of  the  parameters 
fails  to  meet  the  convergence  conditions.  The  difference  is 
that  the  program  will  continue  running  for  this  last  error 
check  and  output  the  results  of  ZSPOW.  The  reason  for  this 
has  already  been  discussed. 

Two  final  comments  on  error  messages  may  help  the 
program  user.  First  of  all,  when  the  convergence  conditions 
are  not  met,  the  program  will  output  an  error  message 
listing  the  specific  case  that  was  not  satisfied.  These 
cases  were  listed  earlier  under  the  subheading  of  phase  two. 
Second,  ZSPOW  outputs  the  following  error  messages  in  Tape  6 

1.  IER=129  indicates  that  the  maximum  number  of 
iterations  has  been  exceeded, 

2.  I£R=130  indicates  that  the  desired  number  of 
significant  digits  is  too  large,  and 

3.  IER=131  indicates  that  ZSPOW  has  not  made  good 
progress . 

All  the  error  checks  were  verified  by  making  the  appropriate 
mistake  on  the  input  tape. 

The  last  comments  added  to  the  program  consisted  of  the 
data  input  format  for  Tape  8.  Although  many  of  the  terms 
have  been  defined  previously  throughout  Chapter  5,  the  defi¬ 
nitions  are  restated  so  they  may  be  seen  in  one  convenient 
location.  If  the  zeroth  moment  is  used,  then  FLAG  =  0.  If 
the  moments  start  with  the  first  moment,  then  FLAG  =  1  .  The 
variable  M  is  the  number  of  "B"  gamma  functions  in  the 


numerator,  N  is  the  number  of  "A”  gamma  functions  in  the 
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numerator,  P  is  the  total  number  of  "A"  gamma  functions,  and 
Q  is  the  total  number  of  "B"  gamma  functions.  If  previously 
calculated  moments  are  used,  then  TYPE  =  0.  If  univariate 
data  are  input,  then  TYPE  =  1.  If  ordered  pairs  from  a 
relative  frequency  or  from  a  continuous  function  are  input, 
then  TYPE  =  2  or  TYPE  =  3  respectively.  GUESS  =  O  if  the 
user  wishes  to  supply  his  own  initial  guess  of  the  H- 
function  parameters.  If  the  default  initial  guess  is  de¬ 
sired,  then  GUESS  =  1.  Finally,  if  TYPE  =  0,  then  NUM  is 
the  number  of  moments.  If  TYPE  =  1,  then  NUM  is  the  number 
of  univariate  data  points.  The  variable  NUM  is  the  number 
of  ordered  pairs  if  TYPE  =  2  or  TYPE  =  3. 

The  data  must  be  input  on  Tape  8  as  follows: 

1.  FLAG,  M,  N,  P,  Q,  TYPE,  GUESS,  NUM  (Integers) 

2.  Data  (Real)  - 

A.  TYPE  0  -  2(P+Q)+2  moments  from  lowest 

to  highest 

B.  TYPE  1  -  NUM  univariate  data 

C.  TYPE  2  or  3  -  NUM  data  pairs  in  the 

form  x,  f(x) 

3.  Initial  Guess  (Real)  - 

A.  GUESS  0  -  2(P+Q)  parameters: 

“B"  pairs  in  the  numerator 
"A"  pairs  in  the  numerator 
"B"  pairs  in  the  denominator 
"A"  pairs  in  the  denominator 

B.  GUESS  1  -  No  input  necessary 


The  program  was  now  in  final  form  and  can  be  seen  in  Appen¬ 


dix  C.  It  had  been  verified  and  validated  by  over  one 
hundred  separate  runsl  We  next  discuss  the  program  limita¬ 
tions  that  came  about  because  of  the  decision  to  use  the 
first  definition  of  the  H-function  in  Chapter  2  and  the 
decision  to  use  the  method  of  moments  in  Chapter  3. 

Program  Limitations 

One  of  the  first  decisions  that  was  made  was  to  use  the 
first  equation  in  Eq  (2.1)  as  our  definition  of  the  H- 
function.  This  decision  forced  the  calculation  of  I(s)  to 
be  done  with  the  formula 

m  n 

TT  rWBiS)  T T  rU-ai-AiS) 

I(s)  =  j=l  3  J  i=l _ 

P  3 

TJ  P(ai+Ais)  TT  rU-b.-B.s) 
i=n+l  j=m+l  J 

This  formula  was  placed  in  a  program  subroutine  named 
COMPIS.  The  program  cannot  solve  for  H-functions  such  as 

m  n 

H ( z )  =  H  [z^ : { ( a^ , A^) } ;  { ( b j , B  j ) } ] 

p  q 


or 


H  ( z ) 


m  n 

H  [I: { (ai,Ai)} ;  {(bj.Bj)}] 

n  n  ^ 


because  the  variable  z  must  have  a  power  of  one.  This  is 
not  a  serious  limitation.  The  program  user  can  use  the 


properties  outlined  in  Chapter  2  in  order  to  convert  the 
H-function  to  the  proper  form  of 


H(z)  =  kH  [cz: { (ai#Ai) } ;  ((bj,Bj)}] 
p  q 


As  an  example,  suppose  the  program  user  expected  an  H- 
function  like 


H(z)  =  H  C(2z)-1:(-l,l);3 
1  0 


He  would  only  be  able  to  run  the  program  after  using  the 
reciprocal  property  to  produce 


H(z)  =  H  £  2z : ; (2 , 1 ) ] 
O  1 


where  the  variable  z  is  taken  to  the  first  power . 

The  second  decision  to  use  the  method  of  moments  posed 
more  serious  problems.  The  first  problem  involved  the 
unstable  nature  of  the  higher  data  moments.  A  sufficient 
discussion  of  this  problem  can  be  found  in  Chapter  3.  With 
data  moments  not  close  to  their  true  value,  the  H-function 
cannot  be  expected  to  accurately  fit  the  true  distribution. 

The  second  problem  with  moments  eliminated  some  func¬ 
tions  and  one  distribution  from  our  H-function  curve-fitting 
procedure.  The  method  of  moments  does  not  apply  to  those 


statistical  distributions  and  functions  that  either  do  not 


have  defined  H-function  moments  (Pr)  or  do  not  have  finite 
data  moments  (Mr).  This  included  all  trigometric  functions, 
the  log  (1-z)  function,  and  the  Half-Cauchy  distribution. 
The  method  of  maximum  likelihood  estimation  may  allow  these 
few  remaining  functions  and  the  one  distribution  to  be  fit 
once  the  theory  has  progressed. 

The  third  problem  involved  restrictions  on  three  of  the 
remaining  2nd  order  distributions.  These  restrictions  come 
from  the  fact  that  2nd  order  H-functions  require  five  addi¬ 
tional  moments  besides  the  zeroth  moment.  Recall  that  the 
recommended  configuration  for  the  program  is  FLAG  =  0.  In 
order  to  generate  the  required  number  of  H-function  moment 
equations,  the  following  restrictions  must  be  met.  For  a 
half-student  distribution,  the  parameter  9  cannot  be  less 
than  or  equal  to  five.  For  an  F  distribution,  the  parameter 
0  cannot  be  less  than  or  equal  to  ten.  Finally,  for  a  beta 
distribution  of  the  second  kind,  the  parameter  0  cannot  be 
less  than  or  equal  to  five. 

The  fourth  and  final  problem  caused  by  the  decision  to 
use  the  method  of  moments  has  already  been  referred  to  as 
the  problem  of  moments.  Moments,  in  general,  do  not 
uniquely  determine  the  distribution  function  when  only  a 
finite  number  of  moments  are  available  (29:81;  33:202-203). 
This  can  be  seen  in  Figure  2,  taken  from  Ramberg  (33:205), 
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1.0 


II 


Figure  2.  Two  Densities  With  Approximately 
the  Same  First  Four  Moments 


which  compares  the  graphs  of  two  different  distributions 
that  have  approximately  the  same  first  four  moments. 

Data  Generation 

The  next  step  was  to  create  data  sets  from  .mathematical 
functions  and  statistical  distributions.  We  decided  not  to 
use  either  type  1  or  type  2  data.  Recall  that  type  1  data 


are  univariate  deviates  and  that  type  2  data  are  ordered 
pairs  from  a  relative  frequency.  Some  mathematical  func¬ 
tions  could  not  be  converted  into  either  type  of  data.  For 
that  matter,  the  conversion  of  univariate  statistical  data 
into  a  relative  frequency  was  both  time-consuminq  and  sub¬ 
jective.  Since  graphs  were  desired,  paired  data  were 


necessary. 


Type  3  data  had  none  of  these  disadvantages.  Recall 
that  type  3  data  are  observations  of  the  function  f(x)  at 
various  values  of  x.  Both  mathematical  functions  and  sta¬ 
tistical  distributions  could  be  easily  converted  into  type  3 
data.  The  data  pairs  made  graphing  possible.  These  graphs 
were  important  to  provide  a  visual  representation  of  the  fit 
of  the  estimated  H-function  to  the  data. 

The  following  three  factors  were  considered  when  we 
created  the  type  3  data: 

1.  Amount  of  data, 

2.  Type  of  function,  and 

3.  Order  of  the  H-function. 

The  amount  of  data  could  either  have  a  small  level  or  a 
large  level.  For  the  small  level,  20  data  points  were  used. 
For  the  large  level,  100  or  150  data  points  were  used.  The 
type  of  function  could  be  either  a  nonstatistical,  mathema¬ 
tical  function  or  a  statistical  distribution.  Examples  of 
mathematical  functions  were  the  power  and  generalized  gamma 
functions.  Statistical  distributions  included  the  exponen¬ 
tial  and  beta  probability  density  functions.  The  third 
factor  also  had  two  levels,  because  we  created  both  first 
and  second  order  data. 

As  a  minimum,  we  created  8  separate  data  sets  to  cover 
each  factor  at  two  levels.  This  demonstrated  the  versatil¬ 
ity  of  the  H-function  curve-fitting  procedure.  In  addition. 


9 


many  other  data  sets  were  created  so  that  more  than  one 
graph  could  be  inspected. 

Procedure 

Once  the  parameter  sets  were  generated,  we  used  the 
computer  program  named  THESIS  (Appendix  C)  to  estimate  the 
parameters  of  the  H-function  from  the  raw  data.  These 
estimates  could  then  be  used  in  two  ways. 

First,  if  the  data  was  from  a  distribution,  then  the 
estimates  could  be  compared  to  Table  III  in  order  to  deter¬ 
mine  the  closest  distribution.  As  an  example,  if  THESIS 
returned  the  values 
b  =  6.00 
B  =  1.01 
k  =  .0027 
c  =  .5 

then  the  program  user  could  hypothesize  that  the  data  was 
from  a  gamma  distribution.  Since  0-l=b  and  b=6.00,  the  user 
would  estimate  9=7.00.  Further,  since  c=0  and  c=.5,  the 
estimate  of  0  is  1/2.  Therefore,  the  user  could  stop  and 
say  the  data  was  from  a  gamma  (0=7,0=.5).  However,  the  user 
implies  B=1.00  and  k=. 0027778  when  he  assumes  that  the  data 
is  from  that  specific  gamma  p.d.f.  There  is  no  general 
method  to  reevaluate  the  other  parameters  once  the 
additional  constraints  have  been  added. 

It  is  better  to  use  the  estimates  of  the  parameters  of 
the  H-function  in  another  way,  as  is  done  in  this  thesis. 
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Once  the  estimated  parameters  were  found,  we  used  another 
computer  program  [7:Appendix  B]  to  determine  the  value  of 


kH  (cx)  for  the  same  values  of  x  as  the  data  pairs. 

p  q 

Finally,  both  groups  of  paired  data  were  graphed  together  to 
display  the  fit  of  the  estimated  H-function  to  the  actual 
data.  These  results,  along  with  the  discussion  of  the 
measure  of  fit  between  the  H-function  and  the  data,  will  be 
described  in  the  next  chapter. 


VI 


Results 


Measure  of  Merit 

Earlier  we  mentioned  the  fact  that  we  would  use  the 
estimated  mean  squared  error  (MSE)  as  our  criterion  for 
measuring  the  goodness-of-fit  of  the  H-function  to  the  data. 
For  the  i*-*1  value  of  x,  let  the  value  for  the  actual  func¬ 
tion  be  f(x^)  and  the  value  for  the  estimated  H-function  be 
H(x^).  Then  the  square  of  the  distance  between  the  H- 
function  and  the  data  point  would  be  [H(x^)-f (x^) ]2.  The 
sum  of  squares  error  (SSE)  is  then  found  by  adding  the 
square  of  the  distance  for  all  n  values  of  x^: 

n 

SSE  =  V  [HUJ-fUi)]2 
i-1  1  1 


Since  we  wanted  to  compare  graphs  with  different  amounts  of 
data,  the  SSE  needed  to  be  adjusted  for  the  number  of  data 
points  as  in 

Estimated  MSE  =  SSE/n 

We  expected  to  see  a  lower  estimated  MSE  for  the  graphs 
generated  from  a  large  amount  of  data. 

Graph  Description 

An  invented  example  of  the  graphs  is  shown  in  Figure  3 
below.  This  graph  was  not  derived  using  the  H-function 
curve- fit  ting  procedure.  As  seen  in  Figure  3,  a  line  will 
be  used  to  represent  the  H-function  data.  The  actual  data 
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Figure  3.  Sample  graph 


will  be  plotted  as  individual  points.  The  estimated  MSE  can 
be  found  in  the  upper  right-hand  corner  of  the  graph. 

The  thesis  procedure  used  input  decks  like  the  follow¬ 
ing  for  a  Gamma  (9  =  3,  0=2)  with  the  amount  of  data  equal  to 


twenty: 


100=0, 

1,0,0,1,3,0,20 

110  = 

.2500 

, 1516326649 

120= 

.3000 

.3678794412 

130= 

.7500 

.5020428603 

140= 

1.0000 

.5413411329 

150= 

1.2500 

.5130312414 

160= 

1.5000 

.4430836153 

170  = 

1.7500 

.3699179469 

180= 

2.0000 

.2930502222 

190= 

2.2500 

.2249571799 

200= 

2.5000 

. 1684486750 

210= 

2.7500 

. 1236243360 

220= 

3.0000 

.0892350784 

230= 

3.2500 

.0635203059 

240  = 

3,5000 

.0446822163 

250  = 

3.7500 

.0311109958 

260= 

4.0000 

.0214696082 

270= 

4.2500 

.0147005397 

280= 

4.5000 

.0099961941 

290= 

4.7500 

, 006755377 6 

300= 

5.0000 

.0045399930 

310=2. 

1234 

320=1 . 

23.45 

Line  100  contains  FLAG,  M,  N,  P,  Q,  TYPE,  GUESS,  and  NUM. 
Lines  110  through  300  contain  the  type  3  data.  The  remain¬ 
ing  two  lines  are  the  initial  guess  of  the  H-function  param¬ 
eters.  Thes’.  data  decks  were  input  into  the  program  named 
THESIS  which  produced 


RESULTS  OF  ZSF'OW 

numerator: 

SMALLER  1  )  = 

B 1 0  B  (  1  >  = 

denominator: 

VALUES  OF  U  C 

K  =  1.112206977428968457 

C=  1 . 102019335254971269 

F  N  0  R M=  .000000000 0 0 0 0 0 0  0 0  0 


results  such  as 

1 . 181637000173736851 
.790076163379332236 

are: 
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FNORM  is  a  measure  of  how  close  the  estimated  root  vector, 
x,  is  to  the  true  solution  of  the  system  of  equations: 


2 ( P+Q ) 

FNORM  =  Y  [F:  (x)  ]‘ 
i=l  1 


where  F^  is  the  it'*1  nonlinear  equation.  These  results  were 
then  converted  into  familiar  H-function  notation  such  as 


(1.11221)  H  C(1.10202)x: ; (1. 18164, .79008)] 
0  1 


Once  this  was  done,  the  estimated  H-function  was  evaluated 
for  the  same  values  of  x  as  the  actual  data.  Then  both  sets 
of  paired  data  were  graphed  in  Figure  4A. 

Other  figures  for  small  and  large  amounts  of  data  were 
created  using  th\  ame  procedure.  They  can  be  found  in  the 
following  lists: 


First  Order  Statistical 

Fig 

ures 

1. 

Gamma  ( 9=3, 0=2) 

4 

a/b 

2. 

Exponential  (0=1/2) 

5 

a/b 

3. 

Chi-Square  (©*=4) 

6 

a/b 

4. 

Weibull  (9=3/2, 0=1) 

7 

a/b 

5. 

Rayleigh  (0=4) 

8 

a/b 

6. 

Maxwell  (©=2) 

9 

a/b 

7. 

Half-Normal  (©=1) 

10 

a/b 

Second  Order  Statistical 


Figures 


1. 

Beta  (9=2, 0=3) 

11 

A/B 

2. 

Power  Function  (9=3) 

12 

A/B 

3. 

Uniform 

13 

A/B 

4. 

Half-Student  (9=16) 

14 

a/b 

5. 

Bessel  (9=1, 0=1) 

15 

First  Order  Functional 

Figures 

1. 

Generalized  Gamma  Function 

( b=l , B=1 / 2 ) 

16 

A/B 

2. 

Generalized  Gamma  Function 

( b=l , B=1 ) 

17 

A/B 

Second  Order  Functional 

Figures 

4.  iu(l-z)Ta  (b=l/2,a=l) 

5.  zb(l-z)+a  (b=l/2,a=2) 


18  . 

19 

20  . 

21  A/B 

22  A/B 


The  figure  number  can  be  used  to  reference  the  correct 
graph.  The  graphs  were  kept  together  and  placed  in  Appendix 
D  for  quick  comparison. 

In  addition  to  Appendix  D,  Table  IV  summarizes  the 
estimated  MSE  for  each  figure.  The  low  estimated  MSE  for 
each  graph  demonstrates  the  ability  of  the  H-funetion  to  fit 
various  sets  of  raw  data.  Another  reassuring  finding  is 
that  the  H-function  fits  better  when  more  data  is  available. 
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TABLE  IV 


Estimated  Mean  Squared  Errors 


Figure 


AMOUNT  OF  DATA 


Small 


Large 


4  A/B 


0.0000658751 


0.0000003710 


5  A/B 


0.0003411134 


0.0000588732 


6  A/B 


7  A/B 


0.0000153797 


0.0000001745 


0.0005039677 


0.0000054808 


8  A/B 


0.0001306474 


0.0000000673 


9  A/B 


0.0000007455 


0.0000000006 


10  A/B 


0.0021639213 


0.0001184473 


11  A/B 


0.0018068877 


0.0000004352 


12  A/B 


0.0027320795 


0.0000056500 


13  A/B 


0.0052858926 


0.0016352035 


14  A/B 


0.0014194577 


0.0002971406 


15 


0.0000537173 


16  A/B 


0.0001038111 


0.0000000003 


17  A/B 


0.0000673331 


0.0000006979 


18  A/B 


0.0002744282 


0.0000006392 


19  A/B 


0.0006307255 


0.0000010364 


20  A/B 


0.6901255958 


2.0283752311 


21  A/B 


0.0029725091 


0.0000013842 


22  A/B 


0.0029537131 


0.0000013089 
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This  can  be  seen  when  the  estimated  MSE  values  for  large 
amounts  of  data  are  lower. 

Demonstration  Runs 

The  list  above  covers  many  statistical  distributions 
and  mathematical  functions.  In  fact,  all  first  order  sta¬ 
tistical  distributions  listed  in  Chapter  2  were  run  (Figures 
4-10).  For  second  order  statistical  distributions,  those 
without  restrictions  on  the  distribution  parameters  were  run 
(Figures  11-13,  15).  In  addition,  the  second  order  half¬ 
student  distribution  was  run  as  a  first  order  distribution 
to  analyze  the  effect  of  underestimating  the  true  order 
(Figure  14).  Finally,  a  group  of  first  and  second  order 
mathematical  functions  with  various  shapes  was  run  (Figures 
16-22). 

All  the  runs  above  were  done  with  type  3  raw  data.  The 
program  could  also  be  used  with  previously  calculated  per¬ 
fect  moments  to  demonstrate  the  special  properties  of  the  H- 
function  discussed  in  Chapter  2.  The  input  decks  differed 
slightly  because  the  type  3  data  was  replaced  by  a  single 
line  of  perfect  moments. 

The  reduction  properties  were  demonstrated  as  follows. 
Perfect  moments  were  derived  for  an  exponential  (0-1/2). 
For  the  first  reduction  method,  an  "A"  gamma  function  in  the 
denominator  had  to  have  the  same  parameters  as  a  "B"  gamma 
function  in  the  numerator.  When  these  two  cancelled,  an 
additional  correct  "B"  gamma  function  had  to  be  left  in  the 
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2  O 

numerator.  The  input  deck  for  the  kH  [cx]  looked  like 

1  2 


1 00=0 t 2 t 0 f 1 t 2 f 0 r 0 • 8 

110=1 ,2, 8, 43, 384, 3840,46080 f 645120 

120=. 342 

130=1.316 

140=2.222 

150=. 106 

160=2.987 

170=. 543 


The  results  from  THESIS  were 


RESULTS  OF  ZSF0W 
NUMERATOR l 
SMALLB  < 1 )  = 
BIGB ( 1 )  = 
SMALLB <  2  >  = 
BIGB(2)= 

denominator: 

SMALLA ( 1 ) = 
BIGA ( 1 ) = 


-.000000107534436133 

.999999863516116250 

2.772153771162535918 

.856711887061141653 

2 . 772153730034659702 
.856711753324468361 


VALUES  OF  K  S  C  AREJ 

K=  .  .499999821510282771 
C=  .499999977412008079 

•  FN0RM=  .000000000000000899 


Note  the  nearly  equivalent  parameters  in  the  numerator  and 
the  denominator.  When  reduced,  the  second  order  H-function 
approximated  the  true  first  order  H-function  of 

IH0  1  (0,1)]. 


Therefore,  the  first  reduction  method  worked 


The  second  reduction  method  works  when  an  or 


vanishes.  For  the  same  exponential  distribution,  the  input 


1  0 

deck  for  kH  [ex]  was 
1  1 


100=0,1,0,1,1,0,0,6 

110=1,2,8,48,384,3340 

120=. 342 

130=1.316 

140=2 . 2'?'? 

150=. 106 


The  results  from  THESIS  were 


RESULTS  OF  ZSP0W  - 
numerator: 

SMALLB  < 1 ) =  .000000000070760199 

BIGB< 1 )=  1.000000000109487530 

denominator: 

SMALL A  < 1 ) =  4.035573643344662287 

BIGA<1>=  .000011372100772754  * 

VALUES  OF  K  i  C  ARE: 

K=  3 . 137655798202459323 

C=  .499992800657203418 

FN0RM=  .000000000000000000 


Note  the  A^  in  the  denominator  is  approaching  zero.  There¬ 
fore,  the  "A"  gamma  function  was  essentially  a  constant 
P(4.036)  term  in  the  denominator.  This  term  was  evaluated 
to  be  6.279.  When  brought  out  of  the  integral  and  under  the 
k,  the  true  k  became  3.138  divided  by  6.279  or  0.4998.  This 
approximated  the  true  first  order  H- function  so  the  second 


reduction  method  also  worked 


Next,  the  generalization  of  the  uniform  distribution 

1  0 

was  demonstrated.  The  input  decks  for  kH  [ex]  were 

1  1 


100=0,1,0,1 , 1 ,0,0,6 

110=1, .5, .33333333, .25, .2, .1 6666667 

120=. 14567 

130=. 96543 

140=1.23456 

150=. 916432 


1 00=0 ,1,0, 1,1,0, 0,o 

110=1, .5, .33333333, .25. .2, .16666667 

120=. 675 

130=2.111 

140=. 337 

150=1.998 


Note  that  only  the  initial  guess  of  the  H-function  parame¬ 
ters  was  changed.  The  first  initial  guess  gave  results  from 


THESIS  as  follows: 


RESULTS  OF  ZSF'GW  - 

numerator: 

SMALLB  < 1 )  = 
BIGR<1)=  1. 

denominator: 

SMALLA ( 1 ) = 

BIGA< 1 )=  1. 


.000001122240934309 

1.008864174078524911 

1.000006384974028606 
1.0088650027922909 66 


VALUES  OF  K  X  C  ARE? 

K=  1.008864845751389794 

C=  .999996837738926558 

FN0RM=  .000000000000015405 


v  v  v  v  rvwr  k-  n%  rv> 


These  results  were  close  to  the  nongeneralized  uniform 


distribution  H  [x:(l(l);  (0,1)]. 

1  1 

When  the  initial  guess  was  changed,  the  results  became 


RESULTS  OF  ZSPOW  - 

numerator: 

SMALLB  <  1 )  = 

BIGB  < 1 >  =  15. 

denominator: 

SMALLA  < 1 >  = 

BIGA< 1 )=  15. 


.000082785039857461 

15.426639569366727756 

1 .000092670887106294 
15.426641384669039780 


VALUES  OF  K  S  C  ARE! 

K=  15.427045040837867873 

C=  .999988515673520429 

FN0RM=  .00000000000000391? 


These  results  demonstrated  the  generalized  uniform  distribu¬ 


tion  uH  [x: (l,u) ;  (0,u)]  where  u*15.427.  Looking  back 
1  1 

at  the  first  uniform,  this  also  was  generalized  but  u=1.009. 

When  these  two  specific  H- functions  were  evaluated  for 
x  values,  the  following  results  were  obtained: 


v.vy.v 


*  .  “  •"  «_■  O  ',  o  r 


IW'JV 


U=1 

.009 

u=15 

.427 

.0500 

,999999 

.0500 

.999999 

.1000 

1.000000 

.1000 

1.000000 

.1500 

1.000000 

.1500 

1 .000000 

.2000 

1.000000 

.2000 

1.000000 

.2500 

1 .000000 

.2500 

1.000000 

.3000 

.999999 

.3000 

1.000000 

.3500 

.999999 

.3500 

.999999 

.4000 

.999999 

.4000 

.999999 

.4500 

.999998 

.4500 

.999993 

.5000 

.999998 

.5000 

.999997 

.5500 

.999998 

.5500 

.999996 

.6000 

.999997 

.6000 

.999995 

.6500 

♦999996 

.6500 

.999994 

.7000 

.999995 

.7000 

.999992 

.7500 

.999995 

♦  7500 

,999991 

.8000 

.999993 

.3000 

.999939 

.8500 

.999991 

.  8500 

.999937 

.9000 

.999989 

.9000 

.999935 

.9300 

.999934 

.9500 

.999982 

1,0000 

.999959 

1.0000 

.999979 

Both  H-functions  were  a  good  approximation  to  the  uniform 
distribution . 

The  generalization  of  the  Pareto  distribution  (0=6)  was 

0  1 

also  demonstrated.  The  input  deck  for  kH  [cx]  was 

1  1 


100*0,0,1,1,1,0,0,6 

110=1,1.2,1.5,2,3,6 

120=-5 . 654 

130=1.112 

140=-7 .174 

150=. 987 


This  gave  results  from  THESIS  of 


RESULTS  OF  ZSF'OW  - 

numerator: 

SMALL*  ( 1 )  *  -3 . 83278965570049479  :l. 

BIGA  < 1 )  =  .690398445680816764 


DENOMINATOR  J 

SMALLB( 1 ) =  -4.832794475434639025 

BIGB  <  1 )  =  .690398684011771735 


VALUES  OF  K  S  C  ARE: 

K=  4.142424973174115621 

C=  1.000001072468101881 

FN0RM=  .000000000000014109 


These  results  demonstrated  the  generalized  Pareto  distribu- 
O  1 

tion  u0H  [xs (l-u( 1+0) ,u) ;  (-u(l+0)fu)]  where  u=0.6904. 

1  1 

The  H-function  was  then  evaluated  for  x  values.  These 
data  pairs  were  then  compared  to  an  actual  Pareto 
distribution  with  0*6: 


1.2500 

1.250287 

1.5000 

.351166 

1.7500 

.119367 

2.0000 

.046875 

2.2500 

.020553 

2 . 5000 

.009830 

2.7500 

.005045 

3.0000 

.002743 

3.2500 

.001567 

3.5000 

.000933 

3.7500 

.000575 

4.0000 

.000366 

4.2500 

♦000240 

4.5000 

.000161 

4.7500 

.000110 

5.0000 

.000077 

5.2500 

.000055 

5.5000 

.000039 

5.7500 

.000029 

6.0000 

.000021 

x  » 2 »j0Q 

1.258291 

1.5000 

.251166 

1 . 7500 

.119367 

2.0000 

.046875 

2.2500 

.020553 

2.5000 

.009830 

2.7500 

> 005045 

3.0000 

.002743 

3.2500 

.001567 

3.5000 

.000933 

3.7500 

.000575 

4.0000 

.000366 

4.2500 

.000240 

4,5000 

.000161 

4.7500 

.000110 

5.0000 

.000077 

5.2500 

.000055 

5.5000 

.000039 

5.7500 

.000029 

6.0000 

.000021 

- 


% 


L,. 


The  estimated  H-function  produced  a  good  fit  to  the  actual 


distribution . 


Now  that  some  properties  have  been  demonstrated,  the 


usefulness  of  the  program  is  apparent.  One  further  run  may 


emphasize  this  point.  While  attempting  to  produce  an  exam¬ 


ple  for  the  reduction  property,  the  following  input  data 


deck  was  used 


1 00=0  ,1,0. 1,1, 0,0, 6 

110=1, 2, 8, 48, 384, 3340 

120=. 342 

130=1.316 

140=2.222 

150=. 886 


The  results  of  THESIS  were 
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numerator: 

SMALLB< 1 )= 

BIGB( 1 )=  2. 

denominator: 

SMALL A  < 1 )  = 

BIGA ( 1 ) =  1. 


- . 000000002439904 1 1 9 
2 . 000000030682542729 


. 499999998705995097 
1.000000030809658824 


VALUES  OF  K  &  C  ARE  J 

K=  1 .772453868441800751 

C=  2.000000041492725700 

FN0RM=  .000000000000000000 


We  were  expecting  to  see  an  A^  near  zero.  Instead  the 


I  0  i 

program  produced  1.772  H  [2x:(±,l);  (0,2)].  But  the  per- 

II  2 

1  1  0 

feet  moments  were  from  ±  H  [ixs ; (0, 1)3 .  At  first  glance, 

2  n  i  2 


the  two  H- functions  did  not  appear  to  be  equivalent. 


*.  »  -  •  •  i 


.v.V.v.  v. 


B. 


However,  using  rule  6.1.18  from  Abramowitz  and  Stegun 
(2:256)  and  the  fact  that  k=Vrr  ,  the  two  H-functions  were 
shown  to  be  equal: 


iT^H  [2x:(I,l);  (0,2)3 


=  1  /  C<2»I  (2x)"s  ds 

2wi  J c  n|+s) 
r  22s P( s)  Hd+s) 

=  _ ! _ 2 _  (2x)“s  ds 

2'r>1  -'c  V2  fli*  rd+s) 

2 

-  ik  j[  *-s 

=  I55T  I  rts)  ‘I*1’3  ds 


s)  (Ax)”3  ds 
* 


1  1  0  i 

±H  Ox:  (0,1)3 

2  0  1  2 


In  a  less  theoretical  approach,  the  H-function  was 
simply  evaluated  for  values  of  x.  The  results  that  follow 
correspond  exactly  to  an  exponential  distribution  with 
0«l/2: 


I 


.7500 

.343645 

1.5000 

.236183 

2.2500 

.162326 

3.0000 

.111565 

3.7500 

♦  076677 

4.5000 

.052700 

5.2500 

.036220 

6.0000 

.024394 

6.7500 

.017109 

7.5000 

.011759 

8.2500 

.008082 

9.0000 

.005554 

9.7500 

♦003818 

10.5000 

.002624 

11.2500 

.001803 

12.0000 

.001239 

12.7500 

.000852 

13.5000 

.000535 

14.2500 

.000402 

15.0000 

.000277 

This  concludes  the  discussion  about  the  results  of  the 
thesis.  More  detailed  conclusions  about  the  efficiency  and 
effectiveness  of  the  H-function  curve-fitting  procedure 
follow  in  the  next  chapter.  Also,  new  findings  are 


highlighted  and  further  studies  are  recommended 


VII  Conclusions  and  Remarks 


Summary 

We  have  developed  a  procedure  to  estimate  the  param¬ 
eters  of  the  H-function  which  gives  the  best  fit  to  a  set  of 
raw  data.  The  procedure  uses  the  method  of  moments  and  can 
be  used  with  both  mathematical  functions  and  continuous 
statistical  distributions  defined  over  positive  x.  Our 
computer  program  will  accept  univariate  data,  data  pairs,  or 
moments  previously  calculated  from  data.  The  user  has  the 
option  of  using  the  zeroth  moment  or  beginning  with  the 
first  moment.  The  user  can  supply  his  own  initial  guess  of 
the  parameters  or  allow  the  program  to  use  a  default  initial 
guess.  The  program  automatically  checks  the  initial  guess 
and  final  estimate  of  the  H-function  parameters  against  the 
convergence  conditions  of  the  H-function.  If  the  program 
should  stop  before  completion,  it  also  has  other  diagnostic 
checks  built  in  which  will  give  the  user  some  indication  of 
the  error  which  caused  the  program  to  abort.  If  no  errors 
are  found,  the  program  will  output  the  parameters  of  the 
fitted  H-function. 

The  method  of  moments  does  not  always  produce  the 
"best"  estimates  of  a  distribution's  parameters.  High 
moments  calculated  from  data  tend  to  be  inaccurate.  Fur¬ 
ther,  moments  do  not  uniquely  define  a  distribution.  Still, 
our  experience  bis  shown  the  method  of  moments  to  be  an 
effective  way  to  estimate  the  parameters  of  the  H-function. 


Because  the  analytic  moments  of  an  H-function  are  easily 
derived  using  the  Mellin  transformation,  the  equations  of 
the  method  of  moments  can  be  simply  written. 


We  used  an  IMSL  routine  named  ZSPOW  to  solve  these 
nonlinear  equations  for  the  unknown  parameters.  ZSPOW  con¬ 
tains  Powell's  quasi-Newton  hybrid  algorithm  for  systems  of 
nonlinear  equations.  This  method  requires  a  reasonably 
close  initial  guess  and  does  not  guarantee  convergence  but 
these  restrictions  are  common  to  most  techniques.  Powell's 
method,  also  provides  super  linear  convergence. 

The  estimated  H-function  parameters  can  be  adjusted 
using  Table  III  if  a  named  statistical  distribution  was 
desired.  Alternatively,  they  could  be  used  as  inputs  to 
another  computer  program  (7:Appendix  B)  which  would  calcu¬ 
late  the  H-function  at  certain  values  of  x  and  plot  the 
probability  density  function  (p.d.f.)  and  cumulative 
distribution  function  (c.d.f.). 

We  tested  our  procedure  using  many  mathematical  func¬ 
tions  and  statistical  distributions.  The  results  were 
impressive  and  are  presented  in  Table  IV  and  Appendix  D. 

Since  many  mathematical  functions  and  statistical  dis¬ 
tributions  are  simultaneously  considered  when  the  H-function 
is  fit  to  a  set  of  data,  fewer  separate  tests  are  required. 
This  generalization  alone  will  increase  the  efficiency  of 
curve-fitting  and  density  estimation. 
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Another  benefit  of  simultaneously  considering  many 
functions  and  distributions  i:  that  there  is  a  higher  proba¬ 
bility  of  finding  the  function  or  distribution  which  "best" 
fits  the  data.  Thus,  effectiveness  increases  when  the 
generalized  method  is  used. 

New  Findings 

We  made  several  contributions  to  the  theory  of  H- 
functions.  We  discovered  a  new  reduction  property  for  H- 
functions.  We  corrected  typographical  errors  in  Mathai  and 
Saxena  (26)  for  arcsin(z),  arctanh(z),  and  log(l±z).  We 
gave  generalized  H-function  formulas  for  the  logarithmic 
function  log(z)  and  power  function  zb.  We  showed  that  the 
Pareto  p.d.f.  can  be  expressed  as  an  H-function.  Finally, 
we  generalized  the  H-function  formulas  for  the  Power  Func¬ 
tion  p.d.f.,  the  Uniform  p.d.f.,  and  the  Pareto  p.d.f. 
These  new  results  were  presented  in  Chapter  2  and  the 
details  that  are  not  obvious  are  presented  in  Appendix  A. 
Recommendations  for  Further  Research 

Because  the  c.d.f.  of  a  random  variable  with  an  H- 
function  distribution  is  simply  one  minus  another  H- 
function,  this  constraint  could  be  added  when  fitting  a 
density.  Fitting  both  the  p.d.f.  and  c.d.f.  could  enhance 
the  power  and  versatility  of  the  H-function  method. 

Further  research  is  also  needed  to  develop  maximum 
likelihood  estimates  for  the  H-function  parameters.  The 
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derivation  was  started  and  can  be  seen  in  Chapter  3.  Our 
results  need  to  be  extended. 

Modifications  to  our  program  could  also  be  made  to 
allow  certain  parameters  to  be  fixed  throughout  the  IMSL 
routine  ZSPOW.  For  example,  the  user  might  find  that  after 
the  first  run  is  made,  B=.97.  He  might  want  to  fix  B=1  and 
solve  for  the  other  parameters,  given  this  constraint.  We 
could  not  provide  this  option  easily  because  ZSPOW  does  not 
allow  a  variable  to  be  altered  within  the  subroutine  that 
defines  the  system  of  nonlinear  equations. 


APPENDIX  A 


Outlines  of  Proofs  of  New  Findings 


The  purpose  of  this  appendix  is  to  outline  proofs  which 
verify  our  new  theoretical  results.  Although  these  proofs 
are  not  central  to  the  purpose  of  the  thesis,  we  felt  com¬ 
pelled  to  verify  our  new  findings  so  that  others  might 
further  extend  the  H- function  theory. 

Throughout  this  appendix  it  is  assumed  that  the  reader 
thoroughly  understands  the  process  of  evaluating  contour 
integrals  in  the  complex  plane  by  summing  the  residues  at 
the  singularities  (or  poles)  of  the  integrand. 

Corrected  Mathematical  Functions 


1  2 


arcsin(z)  *  H  [iz: ( 1 ,1) , ( 1 ,1)  ?  (i,i) ,  (0,1)] 

4f>r  2  2  ‘  2  2  2  2 


E  *  0 


-1  R 


type  VI  +]T  LHP  residues  for  iz)  <  1 
-£  RHP  residues  for  |z|  >  1 

CixiU.i).  (I.i);  (0.1)3 


2  2  '2 
2 


2"  L  rti4.> 


(iz)"s  ds 


numerator  poles  can  be  separated 
by  any  w  in  the  open  interval  (-1,0) 


P(i+is)  has  poles  of  order  1  at 

s j3- ( 2 J ^1 ) i  J =0 f 1| 2f •• • 

(a)  -  (s-sJ)r(Hs)[n“i3)32  (iz)-s 

Ri-is) 

( s+2 J+l )  R  |+|s+ J+l )  [ H  “|s )  3  2  ( iz )  ■ 8 

(i+is)  (i+is+1)  ‘  ’  *  (i+is+J)  ("js)r(-|s) 

» 2n|+^s+j+i)r(“|s)(  iz)-s 

(i+i»)-"(i+is+J-l)(-is) 


_  2f(J+i){iz)2J+1 


(-1)J  Ji  (J+I) 


2iHj+4)z2j+1 

4t» 

Jiu+ij 


arcsin(z) 


_  i  l 
4Y7F  2fTl 


w“2in^>t2Jtl 

J=°  J1  (J+i) 


J  ,4(1.3.5...(2J+l))fSr 

4PW  2J+1  <J+i) 

2 

CD 

I  5]  1.3.5...  ( 2J+1 ) z2J+1 


Z2J+1 


m  z  +  — —  + 


2J  ji  (J+£F 

JL-2  ..Z5  +  — 


2.4.5 


1.3.5  ,7  . 

'»■  -  ■  T  •  •  • 

2. 4. 6. 7 


which  is  the  infinite  series  for  arcsin(z).  The  proof  is 
similar  for  izl  >  1. 


arctanh(z) 


-±  H  [iz«(l,4),  (4,4);  (4,4)  ,  (O, 


2  2 


2  2  2  2 


D  =  1  E  =  0  L  =  -1  R  =  1 
type  VI  +Y,  LHP  residues  for  |z|  <  1 
RHP  residues  for  |z|  >  1 


1  2 


H  [ iz :  ( 1 , 4) ,  (4,4);  (4»i)/  (0,4)3 


2  2 


2  2  2  2 


1/222  2— 2 (iz)“s  ds 


M  L  pu-4s) 


numerator  poles  can  be  separated 

by  any  w  in  the  open  interval  (-1, 

f(4+4s)  has  poles  of  order  1  at 
2  2 

s  ( 2  J  1 ) ,  J **0 , 1 , 2  •  •  • 

2r(4+4a+j+i)nHs)(iz)-s 

fj(3)  58  -T- r1 - r-r— - T — 

(4+4s) . . .  (4+4s+j-*i )  ( -4s) 

2  2  2  2  2 


fjUj)  »  2f~'(  J+l )  (  iz) 2J+1 

(-1)J  Jl  (J+4) 


4i  z 


2J+1 


2J+1 


arctanh( z) 


•  ,  f  S?  „  .  2J+1  ) 

-4  _4_  }  2i>i  T  4l„z  l 

4  2"i  \  jto  2J+1 


0) 


®  2J+1 


M|h- 


which  is  the  infinite  series  for  arctanh(z).  The  proof  is 
similar  for  |z|  >  1. 


log(l±z)  -  H  C±e<(1,1),  (1,1);  (1,1),  (0,1)] 
2  2 

D  »  2  E  »  0  L  *  -1  R  =  1 
type  VI  +£  LHP  residues  for  |z)  <1 
-J]  RHP  residues  for  |z|  >1 


m 


H  [-ZS (1,1),  (1,1);  (1,1),  (0,1)] 
2  2 


f  rd+sjcH-i 

J  c 


UL_  (iz)“s  ds 


numerator  poles  can  be  separated 
by  any  w  in  the  open  interval  (-1,0) 
ru+s)  has  poles  of  order  1  at 
sT=-(J+l),  J=0,l,2... 


£; t(s) 


ndts+Jti)cr(-s)]2(±z)-3 
!i+s) (l+s+i) . . . (l+s+j-iirfi-s! 


\£&A 


ni+stjti)r(-s)(±z)-3 

, l+s)  (1+s+l) ..  .  (1+s+J-lM-s; 


m 


fj(sj) 


P(J-H)(±z)J+1 

(-1)JJI(J+1) 

(±z)J+l 

(-1)J(J+1) 


i"  A-  js'm  u»  ji*'  W|',j»  «■  ■» 


log(l±z)  *  _L_  2-n>i  V  —  t±z)J 

2l>1  J«0  (-1)J(J+1) 

9?  /  XJ+1 

-  y  — <iz> 

J-0  (-1)J(J+1) 

2  3  4 

logU+z)  -  z  -  |  +  |  -  z  +  ... 

2  3  4 

2  3  4 

log ( 1-z)  =  -z  -  f.  -  z  -z  -  ... 

2  3  4 


which  are  the  infinite  series  for  log(l+z)  and  log(l-z) 
The  proof  is  similar  for  jzl  >  l. 

Generalized  Mathematical  Functions 


log ( z )  -  -u2  H  Cz:(l,u),  (l,u);  (0,u),  (0,u)] 
2  2 


O  <  z  <  1 


D®0  E»0  L«-2  R  *  1 


type  VI  +£  LHP  residues  for  |z|  < 
no  RHP  singularities 


log(z) 


z“s  ds 


z~a  ds 


dl 

**  d* 

X  ?  z's  ds 

&  f  -flauL  t-  d. 

**  4  cfd+s)]2 


z”8  ds 
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-H  [z:(l,l),  (1,1);  (0,1),  (0,1)] 

2  2 

0  <  z  <  1 
Q.  E.  D. 


log(z)  “  u2  H  [ z : ( 1 , u ) ,  (l,u);  (0,u),  (0,u)] 

2  2 

D  =  0  E  *  0  L  =  -2  R  =  1 
type  VI  no  LHP  singularities 

-]T  RHP  singularities  for  \z\  >  1 

log(z)  =  — H-r  f  ^-S,n  Z“S  jjg 

2rr:L  Jc  rHl-usH2 


-lii  /  _ 1_  z* 

2tT1  Jc  (-us)2 

i  r  cn-s)]2 

2jrri  Jc  [Rl-s)]2 


c  (-us)2 


z“s  ds 


z“s  ds 


H  Czs(l,l),  (1,1);  (0,1),  (0,1)] 
2  2 


z  >  1 


Q  •  E  •  D  • 


Statistical  Distributions 
Power  Function  p.d.f. 


f(x|e) 


O  <  x  <  1 


0  >  0 


U0H  [x: (u(0-l)+l,u) ;  (u(©-l),u)] 
1  1 

0  E-0  L  *  -1  R  =  1 


U10-1 J+l+US 


2TTi  Jc  (u(©-l)+us)  X  dS 

_®L  f  — * —  x"s  ds 
2-n-i  Jc  0-1+s 


{2iri  x"^1"0)} 

2'fi‘i 


Q.  E.  D. 


Notes  This  generalization  also  applies  to  the  mathematical 
power  function  zb. 

Uniform  p.d.f. 


f(x)  =»  1 


0  <  x  <  1 


1  0 

=  UH  [x: (l,u) ;  (0,u) ] 

1  1 

D*0  E»0  L  =  -1  R  =  1 

type  VI  +Yj  LHP  residues  for  |x|  <  1 
no  RHP  singularities 


«»>  -  *  L  - 

•  “  f  _i  *-■  as 

2TO.  Jc  us 

-ikfj  *-S  - 


I  x“s  ds 


x  ds 


-T-  { 2-n-i  x° }  ■  i 


Q .  E.  D . 


Pareto  p.d.f. 


f(x|9)  =  ex-6-1 


x  >  1 
©  >  0 


O  1 

“  ©H  C  xs (-0, 1 > ;  (-©-1,1)] 

1  1 

D  =  O  E  =  0  L  =  -1  R  =  1 


type  VI  no  LHP  singularities 

-]T  RHP  residues  for  |x|  >  1 


f  ( x  |  ©)  = 


_  © 


2-rri  fc  p(2+©-s) 


jjA±gZ»>  x“s  ds 


.  f  —L- 

2-n-x  L  1+0-s 


x~s  ds 


-  _®_  {2"i  x“(1+e)} 
2«i 


©x 


-0-1 


E.  D. 


Pareto  p.d.f.  (generalized) 

f ( X I  0 )  =  ©X-0-1  X  >  1 

©  >  o 

0  1 

»  u©H  C x: (l-u( 1+©) , u) ;  (-u(l+©),u): 

1  1 

D«0  E  =  0  L  =  -1  R  =  1 
type  VI  no  LHP  singularities 

-Y  RHP  residues  for  |x|  >  1 


APPENDIX  B 


Advanced  Mathematical  Functions  and 
Statistical  Distributions  Expressed  as  H-functions 


The  H-function  also  includes  as  special  cases  several 
advanced  mathematical  functions  and  statistical 
distributions  (26:10-12,145-159?  37:205-206;  7:41-43,87,93). 
Mathematical  Functions 
Bessel  Functions: 


Jv(z) 


4  h1  ° 

z  q  2  *  z  z  z  2 


Kv(z) 


h  q  2  ^  z  z  zz 


Y„(z) 


1  H2  °  l};  (VI)  ,  (-V1),  (-V+11)] 

z  13  Z  ZZ  ZZ  ZZ  2  2 


U  10 

J  (z)  —  H  [z: (0, 1 ) ,  (-V,u)3 
v  0  2 


(Maitland's  generalized  Bessel  function) 


Hypergeoraetric  Functions : 
M(a,b,-z)  =  ^  (a;b;-z) 


1  1 


~  Cr— v  H  [z: (1-a, 1 ) ;  (0,1),  (l-b,l)] 

Pa'  1  2 


(Confluent  Hypergeometric  function) 


2fi  (a,b;c;-z)  =  .J  H  [z:(l-a,l), 

ru)r(o)  2  2 

( 1” b, 1 ) ;  (0,1),  (l-c,l)] 
(Hypergeometric  function) 

rrRbj)  i  p 

pFq  (  ( ai  }  ?  (b  j  } ;  — z)  =  -3Zi -  H  [z  :  {  ( l-ai  ,  1 )  }  ; 

-nTUi)  p  q+1 

i=l 

(0,1),  { ( 1-b j , 1 ) } ] 

for  p<q  or  for  p=q+l  and  I z 1 < 1 . 
(Generalized  Hypergeometric  functions) 


\I/  |"t(ai,Ai)} 

((bi.Bhj*  -z 

*•  J  j  J 


1  P 

=  H  [zsUl-a^Ai)};  (0,1),  { ( 1-b  ■  ,B^ ) }  ] 

p  q+1  J  J 

(Maitland's  or  Wright's  Generalized  Hypergeometric  functions) 


MacRobert's  E- function: 


E(p; {aiJ ;q; {bj} ;z) 

P  1 

=  H  [z:  (1,1),  ((b^Bj);  t  (ajl,Ai) }  3 

q+1  p 


Meijer's  G-function: 


m  n  ,  , 

G  [z: )^i| ] 
P  q  lbj} 


=  H  [z:{(a. ,1)};  ( (b ■ , 1 ) } ] 


Statistical  Distributions 


Bessel  p.d.f. 


f(xls'0>  ~ds  K° 


1 

2rr00 


2  0 
H 

0  2 


(o,i)] 


General  Hypergeometric  p.d.f. 

£  _  c 

f(x|a,b,c,d)  =  dad  T(b)  1  *  _  xc_1M(b,r,-axd) 

r<f>  r<»>-!>  nr) 

=  ad  Rr-d)  H  [adx:  ( l-b+^p»g-)  ; 

rd)r(b-|)  1 2 

(£zi  1),  (i-r+Szi  1)] 


PROGRAM  THESIS 


************************************************************ 


**  ** 

**  WRITTEN  BY:  1LT  RALPH  A.  BOEDIGHEIMER  ** 

**  1LT  CARL  D.  BODENSCHATZ  ** 

**  MS  THESIS  G0R-83D  ** 

**  ** 

**  AIR  FORCE  INSTITUTE  OF  TECHNOLOGY  ** 

**  SCHOOL  OF  ENGINEERING  ** 

#*  WRIGHT-PATTERSON  AIR  FORCE  BASE,  OHIO  ** 

**  ** 

*********  ********************  ******  ************************* 

************************************************************ 
**  ** 

**  PURPOSE:  THIS  PROGRAM  ESTIMATES  THE  PARAMETERS  ** 

**  OF  THE  H-FUNCT I ON  THAT  GIVES  THE  BEST  ** 

**  FIT  TO  A  SET  OF  DATA.  ** 

**  ** 

**  input/output:  tapes/tapeo  ** 

**  ** 

**  NOTE:  THIS  PROGRAM  IS  WRITTEN  IN  FORTRAN  77.  ** 

**  IMSL  MUST  BE  ATTACHED  PRIOR  TO  RUNNING  ** 

**  THE  PROGRAM  FOR  CALLS  TO  ZSPOW  X  GAMMA.  ** 

**  ** 

************************************************************ 

. L»  a,  a.  a.  a.  a>  a.  a.  .a  ti.  a.  a.  a.  a.  a.  tl«  a.  a.  a.  a.  a,  a.  a.  a^  a>  yb  yU  yb  yb  a<  yl#  yk»  ib  ib  a,  dw  a>  yb  vb  «b  ^  vb  yb  yb  yb  yb  yb  a>  yb  yb  yb  a.  yb  a.  ^  yb 

**  ** 

**  FLAG  =  0  IF  THE  ZEROTH  MOMENT  IS  USED  ** 

**  =  1  IF  THE  ZEROTH  MOMENT  IS  NOT  USED  ** 

**  M  =  THE  NUMBER  OF  *B  *  TERMS  IN  THE  NUMERATOR  ** 

**  N  =  THE  NUMBER  OF  "A"  TERMS  IN  THE  NUMERATOR  ** 

**  P  =  THE  TOTAL  NUMBER  OF  *  A  *  TERMS  ** 

**  0  =  THE  TOTAL  NUMBER  OF  *B“  TERMS  ** 

**  TYPE  =  0  IF  MOMENTS  ARE  INPUT  ** 

**  =  1  IF  UNIVARIATE  DATA  ARE  INPUT  ** 

**  =  2  IF  ORDERED  PAIRS  FROM  A  RELATIVE  ** 

**  FREQUENCY  ARE  INPUT  ** 

**  =  3  IF  ORDERED  PAIRS  FROM  A  FUNCTION  ** 

#*  ARE  INPUT  ** 

**  GUESS  =  0  IF  THE  USER  WISHES  TO  SUPPLY  HIS  ** 

**  OWN  INITIAL  PARAMETER  ESTIMATES  ** 

**  =  1  IF  DEFAULT  INITIAL  GUESSES  ARE  DESIRED  ** 

**  NUM  =  THE  NUMBER  OF  MOMENTS  IF  TYPE  =0  ** 

**  =  THE  NUMBER  UF  DATA  POINTS  IF  TYPE  =  1  ** 

**  =  THE  NUMBER  OF  DATA  PAIRS  IF  TYPE  =  2  OR  3  ** 

**  ** 


yb  ^  ^  ^  ^  yb  a.  yb  yb  ^  ^  a.  yb  a*  a.  ^  ^  yb  ^  yb  yb  ^  ^  yb  yb  yb  yb  a.  yb  yb  yb  a.  yb  ^  yb  yb  yb  *b  yb  yb  d.  yb  il»  b  yb  tb  b  yb  yb  yb  yb  ^  yb  yb  ^  yb  yb  yb  ^ 

^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  (p  fp  ^  fp  ^  #p  ^  ^  ^  ^p  ^  ^p  ^  ^  ^p  ^  ^  ^P  ^p  ^  *p  ^p  ^p  ^p  ^  ^  ^P  ^p  ^  ^  ^  ^ 


****  *  *****  *  *  *  *****  :;<*  *  *  **  **  **  *  *  **.' It  **  * . fc  *  **  *  *  *  *  *  * *  *  *  *  *  *  **  ****  *  * 


THE  DATA  MUST  BE  INPUT  IN  TAPE8  IN  THE  FORM: 

1)  FLAG, M,N,P,Q, GUESS, NUM  (ALL  INTEGER) 

2)  DATA  (REAL)  - 

A)  TYPE  0  -  2 ( P+Q ) +2  MOMENTS  FROM 

LOWEST  TO  HIGHEST 

B)  TYPE  1  -  NUM  UNIVARIATE  DATA 

C)  TYPE  2  OR  3  -  NUM  DATA  PAIRS  IN 

THE  FORM  X , F ( X ) 

3)  INITIAL  GUESSES  (REAL)  - 


GUESS 


GUESS  OF 


2 (P+Q)  parameters: 

■B“  PAIRS  IN  NUMERATOR 
•A*  PAIRS  IN  NUMERATOR 
■B’  PAIRS  IN  DENOMINATOR 
*  A  *  PAIRS  IN  DENOMINATOR 
NO  INPUT  NECESSARY 


************************************************************ 


************************************************************ 


EQS  =  THE  NUMBER  OF  EQUATIONS  &  UNKNOWNS 
IER  =  THE  NUMBER  OF  ANY  ERROR  MESSAGE 
ITMAX  =  THE  MAXIMUM  NUMBER  OF  ITERATIONS 
NSIG  *  THE  NUMBER  OF  SIGNIFICANT  DIGITS 
FNORM  =  THE  NORM  OF  THE  F  EQUATION  VECTOR 
PAR  =  A  VECTOR  CONTAINING  FLAG,M,N,P,Q, 
AND  THE  2 ( P+Q ) +2  MOMENTS 
X  =  A  VECTOR  OF  VARIABLES  BEING  ESTIMATED 
(I.E.  THE  ■A*  AND  "B*  PAIRS) 


tl,  tb  ^lj  tb  tb  tb  |b  ab  i b  |b  ab  ab  |b  ab  ^  ab  ala  fb  fb  fb  fb  da  ab  da  ab  a b  \b  a b  (b  a b  fb  ^  da  tb  ajb  |b  ab  ^  ^  |b  df  fb  'b  ^  ^  ^  (b  ^b  df  |b  ^  fb  ^b  da  \b  |b  >b  vb 


INTEGER  CC, EQS, IER, ITMAX, NSIG, FLAG, M,N,P,Q, 

+  TYPE , GUESS , NUM ,I,J,K,L,T,U,V,W 

REAL  FNORM,PAR( 0 t 16 ) ,WK( 418) ,X( 10) , PAIR (4000,2) , 
+  DEV(4000) ,DELTAX, TEST, START1 , START2 , ANSWER 


Uf  ^  df  |b  ^  ^  vb  'b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  ^  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b 

^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  fp  ^  ^  ^  *P  ^  ^  ^  ^  ^  ip  ^  ^  ^  ^  ^  ^  ^  ^  ^  *p  ^  ^  ^  Ip  ^  ^  ip  ^ 

**  ** 

**  WK  IS  A  WORK  VECTOR  WHOSE  SIZE  IS  ** 

**  DEFINED  BY  THE  FOLLOWING  FORMULA:  ** 

**  SIZE  =  ( EQS/2 ) * ( ( 3*EQS ) +15 )  ** 

**  ** 

b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b  b 

^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ®  a^  4a  Ip  ^  ^  ^  1^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ip  ^  ^  ®  ^  1^  ip  4a  ^  ip  ®  ^  ^  Ip  Ip  ®  *  a^  ^ 


EXTERNAL  FCN,COMPIS 


##***##**#**##******##***#****#**  .< * * * * * * * * * * * * * 

%%  %X. 

**  ALL  VARIABLES  ARF  INITIAL I  ZED.  ** 

**  ** 

a************************************************ 

READ  (8,*>  PAR  (  0 )  .PAR  ( 1 )  ,  PAR \ 2) , PAR ( 3 )  , PAR  ( ') )  , 

+  TYPE, GUESS  ,  NUM 

FLAG=PAR ( 0 ) 

M— F’AR  ( 1 ) 

N=PAR ( 2 ) 

P=PAR ( 3 ) 

G-PAR (A) 

EQS=2* ( P+G ) 

NSIG-5 

ITMAX=200 

*****  *  *  *  *  ******  *  *************  *  ****  *  *  **  *  *  *  *  *  *  ****  ** :ic  *  **  * 
**  ** 

**  THE  INPUT  PARAMETERS  OF  FLAG,M,N,P,  AND  G  ** 

**  ARE  CHECKED  AGAINST  SOME  CONDITIONS.  ** 

**  ** 

\L  ilr  tii  Ui  tii  ti#  tli  ti#  til  til  til  til  -ir  til  -ii  tli  tli  tli  t||  tli  tli  tii  -Ji  til  \i/  til  til  tli  til  tii  tii  tli  •  ii  \ii  ti/  U|  u/  -Ji  ii  til  til  vl/  til  til  til  tli  til  *Ji  til  tii  tL  til  tii 

)|(i)(<|il^Yj|(i^i)vl^){(})(i)(ifi<y(i)(iKitt<fl1\ijtifi)(i|vi^i)tfi1(l^lflVki|()flJ|(/fij(lf(i|(i|t/|(^Sif})(i)vif.i|vi)tl|(i|(})ti|(i)(}f.i)t>|(i]t 

IF  < <FLAG.NE.O) .AND. (FLAG. NE.l) )  THEN 
PRINT  *,'  FLAG  MUST  BE  EITHER  0  OR  1 . ' 

GO  TO  999 

ELSEIF  < (M.GT.G) .OR.(M.LT.O) )  THEN 

PRINT  *,'  M  MUST  BE  BETWEEN  0  AND  G,  INCLUSIVE.' 
GO  TO  999 

ELSEIF  ( (N.GT.P) .OR. (N.LT.O) >  THEN 

PRINT  *,'  N  MUST  BE  BETWEEN  0  AND  P,  INCLUSIVE.' 
GO  TO  999 

ELSEIF  ( ( < P+G ) . GT . 5 ) , OR ♦ ( ( P+Q ) , LT , 1 ) )  THEN 
PRINT  *,'  (P+G)  MUST  BE  BETWEEN  1  AND  5.' 

GO  TO  999 
END  IF 


******************************************************** 
**  ** 

**  THE  MOMENTS  ARE  READ  IN  FROM  TAPE  S.  IF  ** 

**  ONLY  RAW  DATA  ARE  AVAILABLE,  THEN  THE  ** 

**  MOMENTS  APE  CALCULATED.  ** 

**  ** 

******************************************************** 


IF  <  TYPE . EQ .0 )  THEN 

IF  <NUM.NE. <EQS+2) >  THEN 

PRINT  THE  NUMBER  OF  MOMENTS  .IS  INCORRECT. 

60  TO  999 
END  IF 

READ  (8,*)  <F*AR<I+4)  ,1  =  1, EQS+2) 

ELSEIF  (TYPE.EQ.l )  THEN 
IF  (NUM.LT .20)  THEN 

PRINT  *,'  FEWER  THAN  20  DATA  POINTS  WILL' 
PRINT  *,'  NOT  PRODUCE  ACCURATE  MOMENTS.' 

GO  TO  999 
ENDIF 

READ  (8,*) <DEV< J) , J=1,NUM> 

CALL  M0M1 <  EOS  ,  FLAG , NUM , DEV , PAR  > 

ELSEIF  < ( TYPE . EQ . 2 ) . OR • ( TYPE  *  EQ . 3 ) >  THEN 
IF  <NUM.LT. 20)  THEN 

PRINT  FEWER  THAN  20  DATA  POINTS  WILL' 

PRINT  NOT  PRODUCE  ACCURATE  MOMENTS.' 

GO  TO  999 
ENDIF 

READ  <3,*>  <  <PAIR<K,L>  ,L.  =  1 ,2)  ,K=1,NUM> 
DELTAX=PAIR< 2,1) -PAIR <1,1 > 

DO  5  J=3,NUM,1 

TEST=PAIR < J , 1 ) -  PAIR < J-l , 1 ) 

IF  <  ABS  <  TEST-BELT  AX ) . GT ♦ < . 5E-5  > )  THEN 
PRINT  DELTA  X  MUST  BE  THE  SAME.' 

GO  TO  999 
ENDIF 
CONTINUE 

IF  <TYPE.EQ.2>  THEN 

CALL  M0M3 < EQS , FLAG , NUM , PAIR , 1 . 0 , PAR ) 

ELSEIF  <TYFE.EQ.3>  THEN 

CALL  M0M3 <  EQS , FLAG , NUM , PAIR , DELTAX , PAR  > 

IF  <FLAG.EQ.O)  THEN 

PRINT  THE  ZEROTH  MOMENT  IS  ',PAR<5> 

PRINT  *,'  IF  YOU  KNOW  THE  DATA  IS  FROM  A' 
PRINT  *,'  STATISTICAL  DISTRIBUTION  AND  WANT 
PRINT  *,'  THE  ZEROTH  MOMENT  TO  EQUAL  ONE,' 
PRINT  *,'  THEN  TYPE  A  ONE.' 

PRINT  *,'  OTHERWISE,  TYPE  ANY  OTHER  NUMBER. 
READ  #, ANSWER 
IF  < ANSWER. EQ. 1.0)  THEN 
PAR<5)=1 .0 
ENDIF 
ENDIF 
ENDIF 
ELSE 

PRINT  #,'  TYPE  MUST  BE  EITHER  0,1,2,  OR  3.' 

GO  TO  999 
ENDIF 


*#  ** 

XX  THE  INITIAL  GUESS  TO  THE  VECTOR  X  IS  READ  IN.  XX 

XX  DEFAULT  VALUES  ARE  AVAILABLE  WHICH  ENSURE  THE  XX 

XX  CONVERGENCE  CONDITIONS  ARE  SATISFIED  INITIALLY  ** 

**  BY  MAKING  D  IN  SUBROUTINE  CHKR  GREATER  THAN  0.  XX 

XX  XX 

XX X X ##*## XX ******** X XXXXXXX X XXXXX * XXXXX X X X X XXX * * * X X X X XXX X X XX 

IF  (GUESS. EQ.O)  THEN 

READ  (8,*> <X(I> ,1=1, EOS) 

ELSEIF  < GUESS ♦ EQ . 1 )  THEN 
START1=. 7654321 
START2=. 87654321 
DO  10  T=1 ,2*M-1 ,2 
X ( T ) =ST  ART1 
X(T+1)=START2 
START1=START1- . 1 
START2=START2- » 1 
10  CONTINUE 

DO  20  U=T,2*(M+N)-1 ,2 
X ( U  >  =ST  ART 1 
X(U  +  .1)=START2 
START 1 =ST ART 1-. 1 
START2=START2-.l 
20  CONTINUE 

DO  30  V=U,2*(Q+N>-1,2 
X ( V ) =ST  ART 1 
X<  V-Fl  )=START2 
ST  ART1=ST  ART1- ♦ 1 
START2=START2- . 1 
30  CONTINUE 

DO  40  U=Vr2*(P+G)-l ,2 
X(W)=START1 
X(U+1)=START2 
START1=3TART1~.1 
START2=ST  ART2- . 1 
40  CONTINUE 

ELSE 

PRINT  VARIABLE  GUESS  MUST  *BE  EITHER  0  OR  1.' 

GO  TO  999 
END  IF 


#  *  *  *  *  *  #  *  *  :K  *  *  #  *  *  *  *  *  *  *  *  *  *  #  *  *  1 1  *  *  *  *  *  *  *  ft  *  *  *  ##  *  *  *  $  :K  *  *  *  * :;:  *  *  #5jc 

**  ** 

**  THE  INITIAL  GUESS  OF  THE  VECTOR  X.  IS  ** 

**  CHECKED  AGAINST  THE  CONVERGENCE  CONDITIONS.  ** 

**  ** 

^  ^  ^  ib  dr  ^  ^  ^  dr  df  ^  ib  lb  d«  ^  ^  ^  ^  lb  ^  lb  lb  d*  ^  ^  vb  ib  ib  |b  d  d/  lb  dr  lb  Os  <b  ib  d*  ib  dr  ib  w  d  *1/  dr  W  *b  lb  ib  vb  ib  d*  ib  ^ 

CALL  CHKR(EQS,M,N,P,Q,X,CC> 

IF  (CC.EQ.O)  THEN 

PRINT  *,'  A  NEW  INITIAL  GUESS  IS  NEEDED.' 

GO  TO  999 
END  IF 

PRINT  *,'  THE  INITIAL  GUESS  MEETS  THE' 

PRINT  *,'  CONVERGENCE  CONDITIONS.' 

#*  ** 

**  ZSPOW  IS  AN  IMSL  ROUTINE  WHICH  USES  ** 

**  POWELL'S  METHOD  TO  APPROXIMATE  THE  ** 

**  ROOTS  OF  A  SYSTEM  OF  EQUATIONS.  *# 

#*  ** 

CALL  ZSPOW  <  FCN , NS I G , EQS , I TMAX , PAR , X , FNORM , WK , I ER ) 

*  *  *###  *  *  #*  #  *#*###  $  *  *  #  *  *  $  #*  xc  #  *##*!((  $  #  *  *****  #**  Xc####*#  Xofc 

**  *# 

**  THE  FINAL  ESTIMATE  OF  THE  VECTOR  X  IS  ** 

**  CHECKED  AGAINST  THE  CONVERGENCE  CONDITIONS.  ** 

*#  ** 


CALL  CHKR(EQS,M,N,P,Q,X,CC) 

IF  (CC.EQ.O)  THEN 

PRINT  *, '  THE  FINAL  ESTIMATE  OF  THE  X  VECTOR  DOES' 
PRINT  *,'  NOT  MEET  THE  CONVERGENCE  CONDITIONS.' 

GO  TO  998 
ENDIF 

PRINT  *,'  THE  FINAL  ESTIMATE  OF  THE  X  VECTOR  MEETS' 
PRINT  *,'  THE  CONVERGENCE  CONDITIONS.' 

d«  >b  df  d.  d/  d.  ib  ib  \b  d/  df  ib  ib  ib  ib  ib  d#  ib  dr  dr  ib  d/  d1  vb  ib  ib  vb  >i.  d/  dr  vb  dr  dr  ib  dr  vb  dr  ib  dr  dr  dr  dr  dr  dr  vb  dr  dr  dr  vlr  dr  dr  dr  dr  ib 

**  ** 

#*  IF  THE  INITIAL  GUESS  OF  THE  X  VECTOR  ** 

**  MEETS  THE  CONVERGENCE  CONDITIONS f  THEN  ** 

**  THE  RESULTS  OF  ZSPOW  ARE  OUTPUT  TO  TAPE  2.  ** 

**  *# 

i***#*######****#####*#***  #****#**#**##**>*#*#  tf#M****** 

998  CALL  RESULT ( EQS , FLAG, M,N,P,Q,X, PAR, FNORM) 

GO  TO  1000 

121 


•  vViL^L*' .'  ■ 


*  *  *  *  *  *  *  &  *  *  *  *  *  *  *  *  *  *  *  .X  *  *  *  *  *  *  *  *  ******  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * 


**  ** 

**  ALL  ERRORS,  EXCEPT  A  FINAL  ESTIMATE  THAT  ** 
**  DOES  NOT  MEET  CONVERGENCE  CONDITIONS,  ARE  ** 
**  ROUTED  TO  THIS  LINE  WHICH  ENDS  THE  ** 
**  EXECUTION  OF  THE  PROGRAM.  ** 
**  ** 


**************************************************** 


999 

PRINT  *," 
PRINT  *," 

THE  PROGRAM  WILL  BE  STOPPED  BECAUSE" 
OF  THIS  ERROR. " 

1000 

CONTINUE 

END 

a********************************************************** 

*********************************************************** 

*********************************************************** 


122 


•  -TV*.  •  .■».  •  »'»  .N  --  \  \  V\ 


SUBROUTINE  FCN (  X  ,  F  » EQS , PAR ) 


**  ** 

#*  THIS  SUBROUTINE  IS  NEEDED  TO  DEFINE  THE  ** 

#*  SYSTEM  OF  NONLINEAR  EQUATIONS*  ** 

**  ** 


J/  U/  Of  U>  Uf  Of  Of  Of  vU  'V  vV  O’  Of  ^  \V  ^  >0/  vl/  O/  Of  vV  \J/  \1/  vl/  vV  0/  SV  Of  Of  ^  Of  Of  \o  >1/  Of  Of  >1/  Of  'if  Of  Of  Of  Of  >1/  vlf  Of  >1/  Or  Of  vV  Of  Of  Of 


INTEGER  EQS ,  I , FLAG ,M»N,P,Q 

REAL  I0FS( 13) »PAK( 0  JEQS+6) , F < EQS ) , X ( EQS ) 

FLAG=PAR ( 0 ) 

M=F'AR<1> 

N=PAR(2) 

P=PAR(3) 

Q*PAR<4) 

CALL  COMP I S ( EQS , FLAG  , M , N , P , Q , X , I OFS ) 

DO  10  1=1, EQS, 1 

F< I)=(PAR< I+4)#PAR< I+6)*I0FS< I+l+FLAG) *#2 > / 

+  <PAR<I+5)**2*I0FS<I+FLAG)#I0FS<I+2+FLAG>)  ~  1.0 

10  CONTINUE 


RETURN 

END 


-*  V  >  V  V  V*  V  V 


TR 


SUBROUT  INF.  RESUL  r  <  EOS  ,  FLAG  ,  M  ,  N  , !-'  ,  G ,  X  ,  PAR  ,  I- NORM  ) 

*#  *>X 

**  THIS  SUBROUTINE  PRINTS  THE  SOLUTION  OF  ** 

**  THE  SYSTEM  OF  NONLINEAR  EQUATIONS.  ** 

**  C  AND  K  ARE  VARIABLES  ESTIMATED  FROM  X.  ** 

##  #* 

INTEGER  EOS , FLAG , M , N , P , Q , T , U , V , W , CNTA  , CNTB 
REAL  IOFS( 13) yX(EQS) ,PAR ( 0 t EQS+6 ) ,FNORM,C,K 

CNTA=0 

CNTB=0 

WRITE ( 2 , 7  <  *  RESULTS  OF  ZSPOW  -■)') 

WRITE ( 2 , 7  <  *  NUMERATOR:*)') 

DO  10  T=1,2*M-1,2 
CNTB=CNTB+1 

WRITE ( 2 7  7  <  *  SMALLB ( * , II , * )*• , F25.18) 7 )CNTB,X<T> 
WRITE <2, 7  <  *  BIGB (  *  ,11,  *>--*,  F25 . 18 )  7  )CNTB,X<T+1 ) 

10  CONTINUE 

DO  20  U=T,2*(M+N)-ly2 
CNTA=CNTA+1 

WRITE (2, 7  < *  SMALLA ( * , I 1 , * ) = * , F25 ♦ 13 ) ' )CNTA,X(U> 
WRITE ( 2 , 7 < "  BIGA( * , 1 1 , * ) * “ , F25 . 18 ) 7 )CNTA,X<U+1) 

20  CONTINUE 

WRITE.  <  2  *  7  <  “  DENOMINATOR:  *  )  7  ) 

DO  30  V=U,2*<G+N)-1,2 
CNTB=CNTB+1 

WRITE <  2 , 7  <  *  SMALLB <  "  ,  1 1 ,  *  >  =  *  , F25 .  13 ) 7 ) CNTB , X ( V ) 
WRITE < 2 r 7 <  *  BIGB< " , II, * )-* ,F25. 18) 7  >CNTB,X(V+1 ) 
30  CONTINUE 

DO  40  W=V,2*(P+Q)-1 ,2 
CNTA=CNTA+1 

WRITE ( 2 , 7 < "  SMALLA ( ’,11,’)=* ,F25.13> 7 )CNTA,X(W) 
WRITE  (2, 7  <*  BIGA(*yIly*)  =  *  ,F25..1.8)  7  )CNTA,X<W+1 ) 

40  CONTINUE 

CALL  COMPIS ( EQS , FLAG , M , N , P , Q , X , IOFS ) 

C= < PAR < 6-FLAG )*J0F5< 3) >/< PAR < 7-FLAG) *IOFS< 2) ) 

K=PAR  <  6-FLAG ) *C##2/I0FS ( 2 ) 

WRITE ( 2 , 7 ( / , "  VALUES  OF  K  %  C  ARE:*)7) 


WRITE ( 2, 7  < * 
WRITE  <  2 , 7 ( * 
WRITE ( 2 , 7 ( * 

RETURN 

END 


K= " , F25  » 1 8 ) 7 ) K 
C=“ ,F25.18) 7 )C 
FN0RM=*  yF25.18,/,/) 7 ) FNORM 


>.v.v. 


screw 


SUH  ROUT  T.NE  COMP  I S  <  EQS  „  FLAG  ,  M  ,  N  ,  P  „  Q  „  X  ,  J  OFS  ) 


**  ** 

**  THIS  SUBROUTINE  COMPUTES  THE  VECTOR  I<S)J  ** 

**  PRODUCTS  AND  QUOTIENTS  OF  GAMMA  FUNCTIONS  FOR  ** 

**  A  GIVEN  VALUE  OF  S.  ** 

#*  ** 


*********************************************************** 

1 NTEGER  EQS , FLAG ,  M ,  N  ,  P  ,  G  ,  R  ,  S ,  T  ,  U ,  V  ,  W 

REAL  X(  EQS )  f  IOFS  ( 13 )» GAMMA  »BNUM  jfBDEN  fANUM  f  ADEN 

DO  5  R-l > 13  *  1 
IOFS  <  R ) =0  ♦  0 
5  CONTINUE 

DO  50  3= 1 +FLAG  >  EQS+2+FLAG » 1 
BNUM-l *0 
BDEN=1  .  0 
ANIJM=1 . 0 
ADEN=1 .0 

DO  10  T»1,2*M-1,2 

BNUM=BNUM*GAMMA<X<T)+S*X<T+1>  ) 

10  CONTINUE 

DO  20  U=T,2*(M+N>-1,2 

ANUM=ANUM*GAMMA< 1-X<U)-S*X<UM. > ) 

20  CONTINUE 

DO  30  V=U,2*(Q+N>-1,2 

BDEN=BDEN*GAMMA< 1-X < V > -S*X < V+ 1 > ) 

30  CONTINUE 

DO  40  W=V,2*(P+Q)-1,2 

ADEN=ADEN*GAMMA<X(W)+S*X(W+1 ) ) 

40  CONTINUE 

IOFS  (  S  )  *  <  BNIJM*ANUM )  /  <  B DEN# ADEN ) 

50  CONTINUE 


RETURN 


SUBROUTINE  CHKR ( EQS , M , N  r P ,Q, X , CC > 


**  ** 

**  THIS  SUBROUTINE  CHECKS  THE  CONVERGENCE  ** 

**  CONDITIONS  FOR  THE  INITIAL  GUESS  AND  THE  *# 

**  FINAL  ESTIMATE  OF  THE  X  VECTOR.  ** 

**  m 


INTEGER  CC , EOS . M , N , P , Q „ T , U , V , U 

REAL  SIJMSBN  ,  SUMBBN  ,  SUMSAN ,  SUMBAN  ,  SUMSBD  .  SUMBBD  , 

+  SUMSAD , SUMBAD , X ( EQS ) , TEST 1 , TESTS . WLOW * WH 1GH 

+  EULQW  ? EWHI GH , D , E , L 

SUMSBN=0 . 0 
SUMBBN=0 . 0 
SUMSAN =0.0 
SUMBAN=0.0 
SUMSBD=0.0 
SUMBBD=0 . 0 
SUMSAB"0.0 
SUMBAD=0  *0 
WLQW=-10000 . 0 
WHIGH=1000Q . 0 

DO  10  T=1,2*M-1,2 
SUMSBN=SUM3BN+X  <  T ) 

SUMBBN=SUMBBNTX  <  T+l  > 

TEST1=-X(T)/X<T+1> 

IF  ( (TEST 1 -WLOW) ,GT. ( ,5E-5) >  THEN 
WLOW-TEST 1 
END  IF 

10  CONTINUE 

DO  20  U=T,2*<M+N>-1  .2 
SUMSAN=SUMSAN+X <  U ) 

SUMBAN=SUMBAN+X ( U+l ) 

TESTS- <  l-X(LJ)  )/X(U  + 1  ) 

IF  <  <  TEST2-WITIGH )  .  LT  .  (-.5E-5)  )  THEN 
WHIGH=TEST2 
END  IF 

20  CONTINUE 

DO  30  V=U,2*<Q+N>-1 ,2 
SUMSBD-SUMSBDTX ( V ) 

SUMBBD=SUMBBD+X ( V+l ) 

30  CONTINUE 

DO  40  W=V,2*(P+Q>-1 ,2 
SUMSAD=SUMSAD+X(W) 

SUMBAD=SUMBAD+X ( U+l ) 

40  CONTINUE 


D-SUMBAN+SUn  B  B  N  -  SUM  B  AD-SI  J  M  B  B  D 
E= <  SUMBAN+SUiiBAD ) - ( SUMBBN+SUMBBD > 

L= < SUMS6N+SUMSBD >  -  <  . 5*Q ) - ( SUMSAN+SUMSAB ) + ( . 5#P ) 


EULOW=E*WLOU 

EWHIGH=E*WHIGH 

IF  ( < WLQW-WHIGH ) .LT . ( - . 5E-5  > )  THEN 
IF  (D.GT.  (  .5F.-5)  )  THEN 
CC=1 

ELSE IF  (B.LT. (-*5E-5) )  THEN 
CC=0 

PRINT  CASE  5 ' 

ELSEIF  < ( ABS <  U ) ♦ LE  .  ( .5E-5) ) . AND . < E .LT . ( -  * 5E-5 ) ) 
IF  < ( L-EWLOW ) . LT ♦  <  - . UE-5 ) )  THEN 
CC=1 

ELSEIF  ( ( L-EWLQW )  .  GT  .  (  . 5E-5 ) )  THEN 

cc=o 

PRINT  CASE  1' 

ELSE 
CC=1 
END  IF 

ELSEIF  ( ( ABS(D) .LE.  < .5E-5) ) .AND. (E.GT. ( .5E-5) ) ) 
IF  (  (L-EUIHIGH)  .LT.  <-.5E-5>  )  THEN 
CC=1 

ELSEIF  ( (L-EWI  IGH) .GT. ( .5E-5) )  THEN 
CC=0 

PRINT  CASE  2 ' 

ELSE 

CC=1 

ENDIF 

ELSEIF  ( (ABS (ID  .LE. ( .5E-5) > .AND. (ABS(E) ,LE*  < .SE 
+  THEN 

IF  <L»LT ♦ (-.5E-5) )  THEN 
CC=1 

ELSEIF  (L.GT, < ,5E-5)  )  THEN 
CC=0 

PRINT  CASE  3' 

ELSE 

CC=0 

PRINT  CASE  4' 

ENDIF 

ENDIF 

ELSE 

CC=0 

PRINT  CASE  6' 

PRINT  NO  OMEGA  IS  POSSIBLE , ' 

ENDIF 


THEN 


THEN 


5)  )  ) 


RETURN 

END 


r * 


,  -  .  v*  .  .• 


SUBROUTINE  M0M.1.  (  EG  3  ,  FLAG  ,  NUM  ,  DE!-' ,  PAR  ) 


**  XX 
**  THIS  SUBROUTINE  GENERATES  THE  MOMENTS  FOR  ** 
**  UP  TO  4000  UNIVARIATE  DATA  POINTS*  ** 
XX  XX 


♦  *  X  XXX  X  X  XXXXX  X  ***  X  X  X  *  X  X  XXX  *  X  :X  XX  X  X  :jc  X  XX  X  X  X  *>!<*  X  *  X  %  %  %  X  XXX  * 

INTEGER  COUNT ,  I ,  J  , EQS , FLAG , NUM 
REAL  SUM , DEV  < 4000 ) , PAR ( 0  X EQ3+A ) 

COUNT =5 

DO  20  I=FLAG,EQS+1+FLAG,1 
SUM=0 . 0 

DO  10  J-l , NUM  , 1 
SUM-SUM+DEV  ( J )  XX .[ 

10  CONTINUE 

SUM=SUM/NUM 
PAR < COUNT )»SUM 
COUNT "COUNT  + 1 
20  CONTINUE 

RETURN 

END 

SUBROUTINE  MOMS (EQS , FLAG , NUM  , PAIR, BELT AX , PAR ) 

*  X  X  X  X  X  X  X  X  X  X  XXXX  *  *  *  *  *  X  X  XX  X  X  X  X  X  X  *  X  X  X  X  XX  X  #  X  X  X  X  X  X  X  *  X  >X  *  X  :jc  X  X 


XX  XX 
XX  THIS  SUBROUTINE  GENERATES  THE  MOMENTS  FOR  ** 
XX  UP  TO  4000  DATA  PAIRS  <X,F(X>>*  XX 
XX  ** 


INTEGER  COUNT , I , J , EQS , FLAG , NUM 

REAL  SUM, PAIR (4000, 2) , DELTAX , PAR ( 0 X EQS+6 > 

COUN  T=5 

DO  20  I=FLAG,EQS+1+FLAG, 1 
SUM=0.0 

DO  10  J=1 , NUM , 1 

SUM=SUM+PAIR< J,1 >#*I*PAIR< J,2) 

10  CONTINUE 

SUM“SUM#DELT  AX 
PAR  ( COUNT  )=SIJM 
COUNT =C0UNT+1 
20  CONTINUE 


RETURN 
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Figure  7A.  Weibull  (0=3/2,  0=1)  n=2Q 


Figure  7B.  Weibull  (9=3/2,  0=1)  n=150 
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Figure  10B.  Half-Norma 
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Figure  14B.  Half-Student  (0=16)  n=150 
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